diff --git a/cime_config/tests.py b/cime_config/tests.py index 19ef3a4477cf..b2371d06db13 100644 --- a/cime_config/tests.py +++ b/cime_config/tests.py @@ -346,7 +346,7 @@ "ERS_Ld5.TL319_oQU240wLI_ais8to30.MPAS_LISIO_JRA1p5.mpaso-ocn_glcshelf", "SMS_P12x2.ne4pg2_oQU480.WCYCL1850NS.allactive-mach_mods", "ERS_Ln9.ne4pg2_ne4pg2.F2010-MMF1.eam-mmf_crmout", - "SMS_Lh4.ne4_ne4.F2010-SCREAMv1.eamxx-output-preset-1", + "SMS_Lh4.ne4_ne4.F2010-SCREAMv1.eamxx-output-preset-1--eamxx-fixer_debug_output", "SMS_Lh4.ne4pg2_ne4pg2.F2010-SCREAMv1.eamxx-output-preset-1--eamxx-prod", ) }, diff --git a/components/eamxx/cime_config/namelist_defaults_eamxx.xml b/components/eamxx/cime_config/namelist_defaults_eamxx.xml index e91e591b1683..765c935794c5 100644 --- a/components/eamxx/cime_config/namelist_defaults_eamxx.xml +++ b/components/eamxx/cime_config/namelist_defaults_eamxx.xml @@ -193,6 +193,10 @@ be lost if SCREAM_HACK_XML is not enabled. moist + false + true + + false 18 diff --git a/components/eamxx/cime_config/testdefs/testmods_dirs/eamxx/fixer_debug_output/shell_commands b/components/eamxx/cime_config/testdefs/testmods_dirs/eamxx/fixer_debug_output/shell_commands new file mode 100644 index 000000000000..2e0f386e2771 --- /dev/null +++ b/components/eamxx/cime_config/testdefs/testmods_dirs/eamxx/fixer_debug_output/shell_commands @@ -0,0 +1,2 @@ +$CIMEROOT/../components/eamxx/scripts/atmchange homme::enable_energy_fixer_debug_info=true -b +./xmlchange POSTRUN_SCRIPT="$CIMEROOT/../components/eamxx/scripts/check-fixer-output" diff --git a/components/eamxx/cmake/tpls/CsmShare.cmake b/components/eamxx/cmake/tpls/CsmShare.cmake index fd682d0ca4e9..d9e667a4e591 100644 --- a/components/eamxx/cmake/tpls/CsmShare.cmake +++ b/components/eamxx/cmake/tpls/CsmShare.cmake @@ -35,6 +35,8 @@ macro (CreateCsmShareTarget) ${SCREAM_BASE_DIR}/../../share/util/shr_orb_mod.F90 ${SCREAM_BASE_DIR}/../../share/util/shr_strconvert_mod.F90 ${SCREAM_BASE_DIR}/../../share/util/shr_sys_mod.F90 + ${SCREAM_BASE_DIR}/../../share/util/shr_reprosum_mod.F90 + ${SCREAM_BASE_DIR}/../../share/util/shr_reprosumx86.c ) # Process genf90 template files. This adds a custom command (and hence target) for each f90 source # that needs to be built from the genf90 template files listed in GENF90_SOURCE. @@ -57,7 +59,8 @@ macro (CreateCsmShareTarget) $<$,$>:CPRGNU> $<$,$>:CPRINTEL> $<$,$>:CPRCRAY> - $<$,$>:CPRCRAY>) + $<$,$>:CPRCRAY> + EAMXX_STANDALONE) if (${CMAKE_SYSTEM} MATCHES "Linux") diff --git a/components/eamxx/scripts/check-fixer-output b/components/eamxx/scripts/check-fixer-output new file mode 100755 index 000000000000..4a63959113fd --- /dev/null +++ b/components/eamxx/scripts/check-fixer-output @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 + +#""" +#aaa +#""" + +import sys, re, glob, pathlib, argparse, gzip + +from utils import run_cmd_no_fail, expect, GoodFormatter + +############################################################################### +def parse_command_line(args, description): +############################################################################### + parser = argparse.ArgumentParser( + usage="""\n{0} [=] ... +OR +{0} --help + +\033[1mEXAMPLES:\033[0m + \033[1;32m# Run hash checker on /my/case/dir \033[0m + > {0} /my/case/dir +""".format(pathlib.Path(args[0]).name), + description=description, + formatter_class=GoodFormatter + ) + + parser.add_argument( + "case_dir", + help="The test case you want to check" + ) + + return parser.parse_args(args[1:]) + +############################################################################### +def readall(fn): +############################################################################### + with open(fn,'r') as f: + txt = f.read() + return txt + +############################################################################### +def greptxt(pattern, txt): +############################################################################### + return re.findall('(?:' + pattern + ').*', txt, flags=re.MULTILINE) + +############################################################################### +def grep(pattern, fn): +############################################################################### + txt = readall(fn) + return greptxt(pattern, txt) + +############################################################################### +def get_log_glob_from_atm_modelio(case_dir): +############################################################################### + filename = case_dir / 'CaseDocs' / 'atm_modelio.nml' + ln = grep('diro = ', filename)[0] + run_dir = pathlib.Path(ln.split()[2].split('"')[1]) + ln = grep('logfile = ', filename)[0] + atm_log_fn = ln.split()[2].split('"')[1] + return str(run_dir / '**' / f'atm.log.*') + +############################################################################### + +# EAMxx:: energy fixer: T tend added to each physics midlevel 0.000222 +# EAMxx:: energy fixer: total energy before fix 32486509719.486938 +# EAMxx:: energy fixer: rel energy error after fix -7.63290383866757e-18 + +def get_fixer_lines(fn,start_from_line): +############################################################################### + fixer_lines = [] + error_vals = [] + + lines = [] + with gzip.open(fn,'rt') as file: + start_line_found = False + for line in file: + if start_line_found: + lines.append(line) + elif start_from_line in line: + start_line_found = True + + i = 0 + while i < len(lines): + line = lines[i] + i = i+1 + # eamxx hash line has the form "eamxx hash> date=YYYY-MM-DD-XXXXX (STRING), naccum=INT + # The INT at the end says how many of the following line contain hashes for this proc-step + if "EAMxx:: energy fixer: rel energy error" in line: + fixer_lines.append(line) + lline = (line.strip('\n')) + errr = float(lline.split()[8]) + error_vals.append(errr) + + return fixer_lines,error_vals + +############################################################################### +def get_model_start_of_step_lines (atm_log): +############################################################################### + lines = [] + with gzip.open(atm_log,'rt') as file: + for line in file: + if "model beg-of-step time" in line: + lines.append(line) + return lines + +############################################################################### +def check_fixer_output(case_dir): +############################################################################### + + ####################### TOLERANCE + TOL = 1e-12 + + case_dir_p = pathlib.Path(case_dir) + expect(case_dir_p.is_dir(), f"{case_dir} is not a dir") + + # Look for the two atm.log files. + glob_pat = get_log_glob_from_atm_modelio(case_dir_p) + atm_fns = glob.glob(glob_pat, recursive=True) + if len(atm_fns) == 0: + print('Could not find atm.log files with glob string {}'.format(glob_pat)) + return False + if len(atm_fns) == 1: + print('Found output file {}'.format(atm_fns[0])) + + start_line = get_model_start_of_step_lines(atm_fns[0])[0] + + print (start_line) + + # Extract hash lines, along with their timestamps, but ignore anything + # before the line $start_line + f = atm_fns[0] + fixer_lines,errs = get_fixer_lines(f,start_line) + print('Array of rel. energy errors after energy fixer:',errs) + errsa = [abs(ele) for ele in errs] + mmax = max(errs) + print('Abs. max. for the rel. energy errors:',mmax) + + if mmax < TOL: + print('SUCCESS') + return True + else: + print('FAIL, abs. max is less than tolerance, which is ', TOL) + return False + + + +############################################################################### +def _main_func(description): +############################################################################### + success = check_fixer_output(**vars(parse_command_line(sys.argv, description))) + sys.exit(0 if success else 1) + +############################################################################### + +if (__name__ == "__main__"): + _main_func(__doc__) diff --git a/components/eamxx/src/control/atmosphere_driver.cpp b/components/eamxx/src/control/atmosphere_driver.cpp index 4b2b114e88ef..04c17c70c71a 100644 --- a/components/eamxx/src/control/atmosphere_driver.cpp +++ b/components/eamxx/src/control/atmosphere_driver.cpp @@ -12,7 +12,7 @@ #include "share/util/eamxx_timing.hpp" #include "share/util/eamxx_utils.hpp" #include "share/io/eamxx_io_utils.hpp" -#include "share/property_checks/mass_and_energy_column_conservation_check.hpp" +#include "share/property_checks/mass_and_energy_conservation_check.hpp" #include "ekat/ekat_assert.hpp" #include "ekat/util/ekat_string_utils.hpp" @@ -420,13 +420,12 @@ void AtmosphereDriver::setup_column_conservation_checks () { // Query m_atm_process_group if any process enables the conservation check, // and if not, return before creating and passing the check. - if (not m_atm_process_group->are_column_conservation_checks_enabled()) { + if (not m_atm_process_group->are_conservation_checks_enabled()) { return; } auto phys_grid = m_grids_manager->get_grid("physics"); const auto phys_grid_name = phys_grid->name(); - // Get fields needed to run the mass and energy conservation checks. Require that // all fields exist. EKAT_REQUIRE_MSG ( @@ -467,7 +466,7 @@ void AtmosphereDriver::setup_column_conservation_checks () const auto heat_flux = m_field_mgr->get_field("heat_flux", phys_grid_name); auto conservation_check = - std::make_shared(phys_grid, + std::make_shared(m_atm_comm,phys_grid, mass_error_tol, energy_error_tol, pseudo_density, ps, phis, horiz_winds, T_mid, qv, diff --git a/components/eamxx/src/dynamics/homme/eamxx_homme_process_interface.cpp b/components/eamxx/src/dynamics/homme/eamxx_homme_process_interface.cpp index dcf48aef7eee..6fc45835b151 100644 --- a/components/eamxx/src/dynamics/homme/eamxx_homme_process_interface.cpp +++ b/components/eamxx/src/dynamics/homme/eamxx_homme_process_interface.cpp @@ -254,7 +254,18 @@ void HommeDynamics::set_grids (const std::shared_ptr grids_m // Create separate remapper for initial_conditions m_ic_remapper = grids_manager->create_remapper(m_cgll_grid,m_dyn_grid); -} + + // Layout for 2D (1d horiz X 1d vertical) variable + const auto& grid_name = m_phys_grid->name(); + // Boundary flux fields for energy and mass conservation checks + if (has_energy_fixer()) { + add_field("vapor_flux", pg_scalar2d, kg/(m2*s), grid_name); + add_field("water_flux", pg_scalar2d, m/s, grid_name); + add_field("ice_flux", pg_scalar2d, m/s, grid_name); + add_field("heat_flux", pg_scalar2d, W/m2, grid_name); + } + +}//set_grids size_t HommeDynamics::requested_buffer_size_in_bytes() const { @@ -468,7 +479,8 @@ void HommeDynamics::initialize_impl (const RunType run_type) // Initialize Rayleigh friction variables rayleigh_friction_init(); -} + +}//initialize_impl void HommeDynamics::run_impl (const double dt) { @@ -747,11 +759,21 @@ void HommeDynamics::homme_post_process (const double dt) { // Store T at end of the dyn timestep (to back out tendencies later) T_prev(ilev) = T_val; }); - }); + }); //op() // Apply Rayleigh friction to update temperature and horiz_winds rayleigh_friction_apply(dt); -} + + if (has_energy_fixer()) { + + get_field_out("vapor_flux").deep_copy(0); + get_field_out("ice_flux").deep_copy(0); + get_field_out("water_flux").deep_copy(0); + get_field_out("heat_flux").deep_copy(0); + + }; //if fixer + +}//homme_post_proc void HommeDynamics:: create_helper_field (const std::string& name, @@ -1048,6 +1070,9 @@ void HommeDynamics::restart_homme_state () { // Erase also qv_prev_phys (if we created it). m_helper_fields.erase("qv_prev_phys"); + + // Update the time stamp of the fields we inited in here (to avoid triggering invalid output in IO) + get_field_out("pseudo_density",pgn).get_header().get_tracking().update_time_stamp(start_of_step_ts()); } void HommeDynamics::initialize_homme_state () { @@ -1223,6 +1248,17 @@ void HommeDynamics::initialize_homme_state () { // Can clean up the IC remapper now. m_ic_remapper = nullptr; + + // Update the time stamp of the fields we inited in here (to avoid triggering invalid output in IO) + get_field_out("pseudo_density",rgn).get_header().get_tracking().update_time_stamp(start_of_step_ts()); + get_internal_field("v_dyn").get_header().get_tracking().update_time_stamp(start_of_step_ts()); + get_internal_field("dp3d_dyn").get_header().get_tracking().update_time_stamp(start_of_step_ts()); + get_internal_field("ps_dyn").get_header().get_tracking().update_time_stamp(start_of_step_ts()); + get_internal_field("phis_dyn").get_header().get_tracking().update_time_stamp(start_of_step_ts()); + get_internal_field("vtheta_dp_dyn").get_header().get_tracking().update_time_stamp(start_of_step_ts()); + if (not params.theta_hydrostatic_mode) { + get_internal_field("w_int_dyn").get_header().get_tracking().update_time_stamp(start_of_step_ts()); + } } // ========================================================================================= void HommeDynamics:: diff --git a/components/eamxx/src/share/CMakeLists.txt b/components/eamxx/src/share/CMakeLists.txt index 8970495f8b05..f40a1420c929 100644 --- a/components/eamxx/src/share/CMakeLists.txt +++ b/components/eamxx/src/share/CMakeLists.txt @@ -34,7 +34,7 @@ set(SHARE_SRC property_checks/property_check.cpp property_checks/field_nan_check.cpp property_checks/field_within_interval_check.cpp - property_checks/mass_and_energy_column_conservation_check.cpp + property_checks/mass_and_energy_conservation_check.cpp util/eamxx_data_interpolation.cpp util/eamxx_fv_phys_rrtmgp_active_gases_workaround.cpp util/eamxx_time_interpolation.cpp @@ -42,6 +42,7 @@ set(SHARE_SRC util/eamxx_timing.cpp util/eamxx_utils.cpp util/eamxx_bfbhash.cpp + util/eamxx_repro_sum_mod.F90 ) # Append ETI sources (I didn't do it above for clarity of reading) @@ -123,6 +124,9 @@ if (EAMXX_ENABLE_EXPERIMENTAL_CODE) endif() add_library(scream_share ${SHARE_SRC}) +set_target_properties(scream_share PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/modules) +target_include_directories (scream_share PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/modules) + target_include_directories(scream_share PUBLIC ${SCREAM_SRC_DIR} ${SCREAM_BIN_DIR}/src @@ -132,6 +136,8 @@ target_include_directories(scream_share PUBLIC # Used in the data interpolation target_link_libraries(scream_share PUBLIC stdc++fs) +target_link_libraries(scream_share PUBLIC csm_share) + if (GPTL_PATH) target_include_directories(scream_share PUBLIC ${GPTL_PATH}) else() diff --git a/components/eamxx/src/share/atm_process/atmosphere_process.cpp b/components/eamxx/src/share/atm_process/atmosphere_process.cpp index f7f615e0347e..9ae79fa84a69 100644 --- a/components/eamxx/src/share/atm_process/atmosphere_process.cpp +++ b/components/eamxx/src/share/atm_process/atmosphere_process.cpp @@ -1,6 +1,6 @@ #include "share/atm_process/atmosphere_process.hpp" #include "share/util/eamxx_timing.hpp" -#include "share/property_checks/mass_and_energy_column_conservation_check.hpp" +#include "share/property_checks/mass_and_energy_conservation_check.hpp" #include "share/field/field_utils.hpp" #ifdef EAMXX_HAS_PYTHON @@ -16,6 +16,9 @@ namespace py = pybind11; #include #include #include +#include +#include + namespace scream { @@ -68,9 +71,28 @@ AtmosphereProcess (const ekat::Comm& comm, const ekat::ParameterList& params) m_repair_log_level = str2LogLevel(m_params.get("repair_log_level","warn")); // Info for mass and energy conservation checks - m_column_conservation_check_data.has_check = + m_conservation_data.has_column_conservation_check = m_params.get("enable_column_conservation_checks", false); + // Energy fixer + m_conservation_data.has_energy_fixer = + m_params.get("enable_energy_fixer", false); + + // Energy fixer + m_conservation_data.has_energy_fixer_debug_info = + m_params.get("enable_energy_fixer_debug_info", false); + + //either energy fixer or column checks, but not both at the same time + EKAT_REQUIRE_MSG ( + !(m_conservation_data.has_energy_fixer && m_conservation_data.has_column_conservation_check), + "Error! In param list " + m_params.name() + " both enable_energy_fixer and" + " enable_column_conservation_check are on, which is not allowed. \n"); + + EKAT_REQUIRE_MSG ( + !((!m_conservation_data.has_energy_fixer) && m_conservation_data.has_energy_fixer_debug_info), + "Error! In param list " + m_params.name() + " enable_energy_fixer is false and" + " enable_energy_fixer_debug_info is true, which is not allowed. \n"); + m_internal_diagnostics_level = m_params.get("internal_diagnostics_level", 0); #ifdef EAMXX_HAS_PYTHON if (m_params.get("py_module_name",std::string(""))!="") { @@ -139,7 +161,16 @@ void AtmosphereProcess::run (const double dt) { m_start_of_step_ts = m_end_of_step_ts; m_end_of_step_ts += dt_sub; - if (has_column_conservation_check()) { + //energy fixer needs both mass and energy fields that are computed by compute_column_conservation_checks_data(dt_sub) + //but this will change with cp* (with cpdry heat_glob const is much easier to compute) + //actually we do not need mass_before, so mass_before is redundant + + //however, we do not want to use column checks if the fixer is on. without a correction + //from the fixer, column checks would fail after fixer. + //we decided that using the column check to verify the fixer (after intorducing a new "flux") + //is not that practical. the fixer will be verified by another call to global integral, + //but disabled by default. + if (has_column_conservation_check() || has_energy_fixer()) { // Column local mass and energy checks requires the total mass and energy // to be computed directly before the atm process is run, as well and store // the correct timestep for the process. @@ -161,6 +192,11 @@ void AtmosphereProcess::run (const double dt) { m_end_of_step_ts, false, true, true); + if (has_energy_fixer()){ + const bool & debug_info = has_energy_fixer_debug_info(); + fix_energy(dt_sub, debug_info ); + } + if (has_column_conservation_check()) { // Run the column local mass and energy conservation checks run_column_conservation_check(); @@ -546,8 +582,8 @@ void AtmosphereProcess::run_column_conservation_check () const { m_atm_logger->debug("[" + this->name() + "] run_column_conservation_check..."); start_timer(m_timer_prefix + this->name() + "::run-column-conservation-checks"); // Conservation check is run as a postcondition check - run_property_check(m_column_conservation_check.second, - m_column_conservation_check.first, + run_property_check(m_conservation.second, + m_conservation.first, PropertyCheckCategory::Postcondition); stop_timer(m_timer_prefix + this->name() + "::run-column-conservation-checks"); m_atm_logger->debug("[" + this->name() + "] run_column-conservation_checks...done!"); @@ -847,11 +883,11 @@ add_postcondition_check (const prop_check_ptr& pc, const CheckFailHandling cfh) void AtmosphereProcess:: add_column_conservation_check(const prop_check_ptr &prop_check, const CheckFailHandling cfh) { - EKAT_REQUIRE_MSG(m_column_conservation_check.second == nullptr, + EKAT_REQUIRE_MSG(m_conservation.second == nullptr, "Error! Conservation check for process \""+ name() + "\" has already been added."); - m_column_conservation_check = std::make_pair(cfh,prop_check); + m_conservation = std::make_pair(cfh,prop_check); } void AtmosphereProcess::set_fields_and_groups_pointers () { @@ -1183,20 +1219,45 @@ ::remove_group (const std::string& group_name, const std::string& grid_name) { rmg(m_groups_out, m_groups_out_pointers); } -void AtmosphereProcess::compute_column_conservation_checks_data (const int dt) +void AtmosphereProcess::compute_column_conservation_checks_data (const double dt) { - EKAT_REQUIRE_MSG(m_column_conservation_check.second != nullptr, + EKAT_REQUIRE_MSG(m_conservation.second != nullptr, "Error! User set enable_column_conservation_checks=true, " - "but no conservation check exists.\n"); + "or has_energy_fixer=true, " + "but no conservation check class exists.\n"); // Set dt and compute current mass and energy. const auto& conservation_check = - std::dynamic_pointer_cast(m_column_conservation_check.second); + std::dynamic_pointer_cast(m_conservation.second); conservation_check->set_dt(dt); conservation_check->compute_current_mass(); conservation_check->compute_current_energy(); } +void AtmosphereProcess::fix_energy (const double dt, const bool &print_debug_info) +{ + EKAT_REQUIRE_MSG(m_conservation.second != nullptr, + "Error! User set has_energy_fixer=true, " + "but no conservation check class exists.\n"); + + // Set dt and compute current mass and energy. + const auto& conservation_check = + std::dynamic_pointer_cast(m_conservation.second); + + //dt is needed to convert flux to change + conservation_check->set_dt(dt); + conservation_check->global_fixer(print_debug_info); + + if(print_debug_info){ + //print everything about the fixer only in debug mode + m_atm_logger->info("EAMxx:: energy fixer: T tend added to each physics midlevel " + std::to_string( conservation_check->get_pb_fixer() ) ); + m_atm_logger->info("EAMxx:: energy fixer: total energy before fix " + std::to_string( conservation_check->get_total_energy_before() ) ); + std::stringstream ss; + ss << "EAMxx:: energy fixer: rel energy error after fix " << std::setprecision(15) << conservation_check->get_echeck() << "\n"; + m_atm_logger->info(ss.str()); + } +} + void AtmosphereProcess::add_py_fields (const Field& f) { #ifdef EAMXX_HAS_PYTHON diff --git a/components/eamxx/src/share/atm_process/atmosphere_process.hpp b/components/eamxx/src/share/atm_process/atmosphere_process.hpp index 937585f184b0..0631f690cb01 100644 --- a/components/eamxx/src/share/atm_process/atmosphere_process.hpp +++ b/components/eamxx/src/share/atm_process/atmosphere_process.hpp @@ -166,8 +166,8 @@ class AtmosphereProcess : public ekat::enable_shared_from_this - get_column_conservation_check() { - return m_column_conservation_check; + get_conservation() { + return m_conservation; } @@ -271,7 +271,9 @@ class AtmosphereProcess : public ekat::enable_shared_from_this& get_restart_extra_data () { return m_restart_extra_data; } // Boolean that dictates whether or not the conservation checks are run for this process - bool has_column_conservation_check () { return m_column_conservation_check_data.has_check; } + bool has_column_conservation_check () { return m_conservation_data.has_column_conservation_check; } + bool has_energy_fixer () { return m_conservation_data.has_energy_fixer; } + bool has_energy_fixer_debug_info () { return m_conservation_data.has_energy_fixer_debug_info; } // Print a global hash of internal fields (useful for debugging non-bfbness) // Note: (mem, nmem) describe an arbitrary device array. If mem!=nullptr, @@ -536,7 +538,9 @@ class AtmosphereProcess : public ekat::enable_shared_from_this> m_postcondition_checks; // Column local mass and energy conservation check - std::pair m_column_conservation_check; + std::pair m_conservation; // Store data related to this processes conservation check. - struct ColumnConservationCheckData { + struct ConservationData { // Boolean which dictates whether or not this process // contains the mass and energy conservation checks. - bool has_check; + bool has_column_conservation_check; // Tolerance used for the conservation check + // mass or energy or both? rename Real tolerance; + bool has_energy_fixer; + bool has_energy_fixer_debug_info; }; - ColumnConservationCheckData m_column_conservation_check_data; + ConservationData m_conservation_data; // This process's copy of the timestamps (current, as well as beg/end of step) TimeStamp m_start_of_step_ts; diff --git a/components/eamxx/src/share/atm_process/atmosphere_process_group.cpp b/components/eamxx/src/share/atm_process/atmosphere_process_group.cpp index abc9f0d93753..33791037d6ef 100644 --- a/components/eamxx/src/share/atm_process/atmosphere_process_group.cpp +++ b/components/eamxx/src/share/atm_process/atmosphere_process_group.cpp @@ -221,7 +221,7 @@ gather_internal_fields () { } bool AtmosphereProcessGroup:: -are_column_conservation_checks_enabled () const +are_conservation_checks_enabled () const { // Loop through processes and return true if an instance is found. for (auto atm_proc : m_atm_processes) { @@ -230,7 +230,7 @@ are_column_conservation_checks_enabled () const // else continue to the next process. auto atm_proc_group = std::dynamic_pointer_cast(atm_proc); if (atm_proc_group) { - if (atm_proc_group->are_column_conservation_checks_enabled()) { + if (atm_proc_group->are_conservation_checks_enabled()) { return true; } else { continue; @@ -239,7 +239,7 @@ are_column_conservation_checks_enabled () const // If this process is not a group, query enable_column_conservation_checks // and return true if true. - if (atm_proc->has_column_conservation_check()) { + if (atm_proc->has_column_conservation_check() || atm_proc->has_energy_fixer()) { return true; } } @@ -249,7 +249,7 @@ are_column_conservation_checks_enabled () const } void AtmosphereProcessGroup:: -setup_column_conservation_checks (const std::shared_ptr& conservation_check, +setup_column_conservation_checks (const std::shared_ptr& conservation_check, const CheckFailHandling fail_handling_type) const { // Loop over atm processes and add mass and energy checker where relevant @@ -265,38 +265,28 @@ setup_column_conservation_checks (const std::shared_ptrhas_column_conservation_check(), "Error! The ATM process group \"" + atm_proc_group->name() + "\" attempted to enable " - "conservation checks. Should have enable_column_conservation_checks=false for all " - "process groups.\n"); - + "conservation checks. A process group cannot have enable_column_conservation_checks=true. \n"); atm_proc_group->setup_column_conservation_checks(conservation_check, fail_handling_type); continue; } // For individual processes, first query if the checks are enabled. // If not, continue to the next process. - if (not atm_proc->has_column_conservation_check()) { + if (not (atm_proc->has_column_conservation_check() || atm_proc->has_energy_fixer()) ) { continue; } - // Since the checker is column local, require that an atm - // process that enables the check is a Physics process. - EKAT_REQUIRE_MSG(atm_proc->type() == AtmosphereProcessType::Physics, - "Error! enable_column_conservation_checks=true " - "for non-physics process \"" + atm_proc->name() + "\". " - "This check is column local and therefore can only be run " - "on physics processes.\n"); - // Query the computed fields for this atm process and see if either the mass or energy computation // might be changed after the process has run. If no field used in the mass or energy calculate // is updated by this process, there is no need to run the check. const std::string phys_grid_name = conservation_check->get_grid()->name(); - const bool updates_static_energy = atm_proc->has_computed_field("T_mid", phys_grid_name); + const bool updates_internal_energy = atm_proc->has_computed_field("T_mid", phys_grid_name); const bool updates_kinetic_energy = atm_proc->has_computed_field("horiz_winds", phys_grid_name); const bool updates_water_vapor = atm_proc->has_computed_field("qv", phys_grid_name); const bool updates_water_liquid = atm_proc->has_computed_field("qc", phys_grid_name) || atm_proc->has_computed_field("qr", phys_grid_name); const bool updates_water_ice = atm_proc->has_computed_field("qi", phys_grid_name); - const bool mass_or_energy_is_updated = updates_static_energy || updates_kinetic_energy || + const bool mass_or_energy_is_updated = updates_internal_energy || updates_kinetic_energy || updates_water_vapor || updates_water_liquid || updates_water_ice; EKAT_REQUIRE_MSG(mass_or_energy_is_updated, "Error! enable_column_conservation_checks=true for " @@ -319,7 +309,8 @@ setup_column_conservation_checks (const std::shared_ptradd_column_conservation_check(conservation_check, fail_handling_type); - } + }// for (auto atm_proc : m_atm_processes) + } void AtmosphereProcessGroup::add_postcondition_nan_checks () const { @@ -358,7 +349,7 @@ void AtmosphereProcessGroup::add_additional_data_fields_to_property_checks (cons prop_check.second->set_additional_data_field(data_field); } if (proc->has_column_conservation_check()) { - proc->get_column_conservation_check().second->set_additional_data_field(data_field); + proc->get_conservation().second->set_additional_data_field(data_field); } } } diff --git a/components/eamxx/src/share/atm_process/atmosphere_process_group.hpp b/components/eamxx/src/share/atm_process/atmosphere_process_group.hpp index 3d7dfab35252..cc099865d89e 100644 --- a/components/eamxx/src/share/atm_process/atmosphere_process_group.hpp +++ b/components/eamxx/src/share/atm_process/atmosphere_process_group.hpp @@ -2,7 +2,7 @@ #define SCREAM_ATMOSPHERE_PROCESS_GROUP_HPP #include "share/atm_process/atmosphere_process.hpp" -#include "share/property_checks/mass_and_energy_column_conservation_check.hpp" +#include "share/property_checks/mass_and_energy_conservation_check.hpp" #include "control/surface_coupling_utils.hpp" #include "ekat/ekat_parameter_list.hpp" @@ -89,12 +89,12 @@ class AtmosphereProcessGroup : public AtmosphereProcess // Returns true if any internal processes enables // the mass and energy conservation checks. - bool are_column_conservation_checks_enabled () const; + bool are_conservation_checks_enabled () const; // Adds the mass and energy conservation // checks to appropriate physics processes. void setup_column_conservation_checks ( - const std::shared_ptr& conservation_check, + const std::shared_ptr& conservation_check, const CheckFailHandling fail_handling_type) const; // Add nan checks after each non-group process, for each computed field. diff --git a/components/eamxx/src/share/property_checks/mass_and_energy_column_conservation_check.cpp b/components/eamxx/src/share/property_checks/mass_and_energy_conservation_check.cpp similarity index 59% rename from components/eamxx/src/share/property_checks/mass_and_energy_column_conservation_check.cpp rename to components/eamxx/src/share/property_checks/mass_and_energy_conservation_check.cpp index 4821f37024fb..6d78bff660b1 100644 --- a/components/eamxx/src/share/property_checks/mass_and_energy_column_conservation_check.cpp +++ b/components/eamxx/src/share/property_checks/mass_and_energy_conservation_check.cpp @@ -1,4 +1,4 @@ -#include "share/property_checks/mass_and_energy_column_conservation_check.hpp" +#include "share/property_checks/mass_and_energy_conservation_check.hpp" #include "physics/share/physics_constants.hpp" #include "share/field/field_utils.hpp" @@ -7,24 +7,26 @@ namespace scream { -MassAndEnergyColumnConservationCheck:: -MassAndEnergyColumnConservationCheck (const std::shared_ptr& grid, - const Real mass_error_tolerance, - const Real energy_error_tolerance, - const Field& pseudo_density, - const Field& ps, - const Field& phis, - const Field& horiz_winds, - const Field& T_mid, - const Field& qv, - const Field& qc, - const Field& qr, - const Field& qi, - const Field& vapor_flux, - const Field& water_flux, - const Field& ice_flux, - const Field& heat_flux) - : m_grid (grid) +MassAndEnergyConservationCheck:: +MassAndEnergyConservationCheck (const ekat::Comm& comm, + const std::shared_ptr& grid, + const Real mass_error_tolerance, + const Real energy_error_tolerance, + const Field& pseudo_density, + const Field& ps, + const Field& phis, + const Field& horiz_winds, + const Field& T_mid, + const Field& qv, + const Field& qc, + const Field& qr, + const Field& qi, + const Field& vapor_flux, + const Field& water_flux, + const Field& ice_flux, + const Field& heat_flux) + : m_comm (comm) + , m_grid (grid) , m_dt (std::nan("")) , m_mass_tol (mass_error_tolerance) , m_energy_tol (energy_error_tolerance) @@ -35,6 +37,8 @@ MassAndEnergyColumnConservationCheck (const std::shared_ptr& m_current_mass = view_1d ("current_total_water", m_num_cols); m_current_energy = view_1d ("current_total_energy", m_num_cols); + m_energy_change = view_1d ("energy change fixer", m_num_cols); + m_fields["pseudo_density"] = pseudo_density; m_fields["ps"] = ps; m_fields["phis"] = phis; @@ -48,9 +52,22 @@ MassAndEnergyColumnConservationCheck (const std::shared_ptr& m_fields["water_flux"] = water_flux; m_fields["ice_flux"] = ice_flux; m_fields["heat_flux"] = heat_flux; + + //allocate Fields for fixer reductions + using namespace ekat::units; + using namespace ShortFieldTagsNames; + + //this does not exist, why? cause it can be a LEV or a COL array? + //FieldLayout scalar1d = m_grid->get_1d_scalar_layout(); + FieldLayout scalar1d{{COL}, {m_num_cols}}; + FieldIdentifier s1_fid("s1", scalar1d, kg / kg, m_grid->name()); + field_version_s1 = Field(s1_fid); + + field_version_s1.allocate_view(); + } -void MassAndEnergyColumnConservationCheck::compute_current_mass () +void MassAndEnergyConservationCheck::compute_current_mass () { auto mass = m_current_mass; const auto ncols = m_num_cols; @@ -76,7 +93,7 @@ void MassAndEnergyColumnConservationCheck::compute_current_mass () }); } -void MassAndEnergyColumnConservationCheck::compute_current_energy () +void MassAndEnergyConservationCheck::compute_current_energy () { auto energy = m_current_energy; const auto ncols = m_num_cols; @@ -107,7 +124,7 @@ void MassAndEnergyColumnConservationCheck::compute_current_energy () }); } -PropertyCheck::ResultAndMsg MassAndEnergyColumnConservationCheck::check() const +PropertyCheck::ResultAndMsg MassAndEnergyConservationCheck::check() const { auto mass = m_current_mass; auto energy = m_current_energy; @@ -277,9 +294,162 @@ PropertyCheck::ResultAndMsg MassAndEnergyColumnConservationCheck::check() const return res_and_msg; } +void MassAndEnergyConservationCheck::global_fixer(const bool & print_debug_info) +{ + using TPF = ekat::TeamPolicyFactory; + const auto ncols = m_num_cols; + const auto nlevs = m_num_levs; + + //keep dt for fixers with fluxes. + //we do not plan to use fluxes in dycore fixer, but the code is already there + EKAT_REQUIRE_MSG(!std::isnan(m_dt), "Error! Timestep dt must be set in MassAndEnergyConservationCheck " + "before running check()."); + auto dt = m_dt; + + const auto pseudo_density = m_fields.at("pseudo_density").get_view (); + const auto T_mid = m_fields.at("T_mid" ).get_view< Real**> (); + const auto horiz_winds = m_fields.at("horiz_winds" ).get_view(); + const auto qv = m_fields.at("qv" ).get_view (); + const auto qc = m_fields.at("qc" ).get_view (); + const auto qi = m_fields.at("qi" ).get_view (); + const auto qr = m_fields.at("qr" ).get_view (); + const auto ps = m_fields.at("ps" ).get_view (); + const auto phis = m_fields.at("phis" ).get_view (); + + const auto vapor_flux = m_fields.at("vapor_flux").get_view(); + const auto water_flux = m_fields.at("water_flux").get_view(); + const auto ice_flux = m_fields.at("ice_flux" ).get_view(); + const auto heat_flux = m_fields.at("heat_flux" ).get_view(); + + const auto policy = TPF::get_default_team_policy(ncols, nlevs); + + using namespace ekat::units; + using namespace ShortFieldTagsNames; + + auto field_view_s1 = field_version_s1.get_view(); + + auto area = m_grid->get_geometry_data("area").clone(); + auto area_view = area.get_view(); + + //reprosum vars + const Real* send; + Real recv; + Int nlocal = ncols; + Int ncount = 1; + +//unite 1st and 2nd ||4 and repro calls + Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const KT::MemberType& team) { + const int i = team.league_rank(); + const auto pseudo_density_i = ekat::subview(pseudo_density, i); + // Calculate total gas mass (sum dp, no water loading) + field_view_s1(i) = compute_gas_mass_on_column(team, nlevs, pseudo_density_i) * area_view(i); + }); + Kokkos::fence(); + + auto s1_host = field_version_s1.get_view(); + send = s1_host.data(); + + field_version_s1.sync_to_host(); + eamxx_repro_sum(send, &recv, nlocal, ncount, MPI_Comm_c2f(m_comm.mpi_comm())); + field_version_s1.sync_to_dev(); + + total_gas_mass_after = recv; + +//this ||4 needs to be 2 for-loops, with one summing each 4 cols first, serially +//for pg2 grids (if on np4 grids, it would require much more work?) + auto energy_change = m_energy_change; + auto current_energy = m_current_energy; + Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const KT::MemberType& team) { + + const int i = team.league_rank(); + const auto pseudo_density_i = ekat::subview(pseudo_density, i); + const auto T_mid_i = ekat::subview(T_mid, i); + const auto horiz_winds_i = ekat::subview(horiz_winds, i); + const auto qv_i = ekat::subview(qv, i); + const auto qc_i = ekat::subview(qc, i); + const auto qr_i = ekat::subview(qr, i); + const auto qi_i = ekat::subview(qi, i); + + // Calculate total energy + const auto new_energy_for_fixer = compute_total_energy_on_column(team, nlevs, pseudo_density_i, T_mid_i, horiz_winds_i, + qv_i, qc_i, qr_i, ps(i), phis(i)); + Kokkos::single(Kokkos::PerTeam(team),[&]() { + energy_change(i) = compute_energy_boundary_flux_on_column(vapor_flux(i), water_flux(i), ice_flux(i), heat_flux(i))*dt; + field_view_s1(i) = (current_energy(i)-new_energy_for_fixer-energy_change(i)) * area_view(i); + }); + }); + Kokkos::fence(); + + field_version_s1.sync_to_host(); + eamxx_repro_sum(send, &recv, nlocal, ncount, MPI_Comm_c2f(m_comm.mpi_comm())); + field_version_s1.sync_to_dev(); + pb_fixer = recv; + + if(print_debug_info) { + //total energy needed for relative error + Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const KT::MemberType& team) { + const int i = team.league_rank(); + Kokkos::single(Kokkos::PerTeam(team),[&]() { + field_view_s1(i) = current_energy(i) * area_view(i); + }); + }); + Kokkos::fence(); + + field_version_s1.sync_to_host(); + eamxx_repro_sum(send, &recv, nlocal, ncount, MPI_Comm_c2f(m_comm.mpi_comm())); + field_version_s1.sync_to_dev(); + + total_energy_before = recv; + } + + using PC = scream::physics::Constants; + const Real cpdry = PC::Cpair; + pb_fixer /= (cpdry*total_gas_mass_after); // T change due to fixer + + //add the fixer to temperature + Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const KT::MemberType& team) { + const int i = team.league_rank(); + const auto T_mid_i = ekat::subview(T_mid, i); + Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlevs), [&] (const int k){ + T_mid_i(k) += pb_fixer; + }); + });//adding fix to T + + if(print_debug_info){ + Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const KT::MemberType& team) { + + const int i = team.league_rank(); + const auto pseudo_density_i = ekat::subview(pseudo_density, i); + const auto T_mid_i = ekat::subview(T_mid, i); + const auto horiz_winds_i = ekat::subview(horiz_winds, i); + const auto qv_i = ekat::subview(qv, i); + const auto qc_i = ekat::subview(qc, i); + const auto qr_i = ekat::subview(qr, i); + const auto qi_i = ekat::subview(qi, i); + + // Calculate total energy + const auto new_energy_for_fixer = compute_total_energy_on_column(team, nlevs, pseudo_density_i, T_mid_i, horiz_winds_i, + qv_i, qc_i, qr_i, ps(i), phis(i)); + //overwrite the "new" fields with relative change + + Kokkos::single(Kokkos::PerTeam(team),[&]() { + field_view_s1(i) = (current_energy(i)-new_energy_for_fixer-energy_change(i))*area_view(i); + }); + }); + Kokkos::fence(); + + field_version_s1.sync_to_host(); + eamxx_repro_sum(send, &recv, nlocal, ncount, MPI_Comm_c2f(m_comm.mpi_comm())); + field_version_s1.sync_to_dev(); + + echeck = recv/total_energy_before; + } + +};//global_fixer + KOKKOS_INLINE_FUNCTION -Real MassAndEnergyColumnConservationCheck:: +Real MassAndEnergyConservationCheck:: compute_total_mass_on_column (const KT::MemberType& team, const int nlevs, const uview_1d& pseudo_density, @@ -302,7 +472,23 @@ compute_total_mass_on_column (const KT::MemberType& team, } KOKKOS_INLINE_FUNCTION -Real MassAndEnergyColumnConservationCheck:: +Real MassAndEnergyConservationCheck:: +compute_gas_mass_on_column (const KT::MemberType& team, + const int nlevs, + const uview_1d& pseudo_density) +{ + using RU = ekat::ReductionUtils; + using PC = scream::physics::Constants; + const Real gravit = PC::gravit; + + return RU::parallel_reduce(team, 0, nlevs, + [&] (const int lev, Real& local_mass) { + local_mass += pseudo_density(lev)/gravit; + }); +} + +KOKKOS_INLINE_FUNCTION +Real MassAndEnergyConservationCheck:: compute_mass_boundary_flux_on_column (const Real vapor_flux, const Real water_flux) { @@ -313,7 +499,7 @@ compute_mass_boundary_flux_on_column (const Real vapor_flux, } KOKKOS_INLINE_FUNCTION -Real MassAndEnergyColumnConservationCheck:: +Real MassAndEnergyConservationCheck:: compute_total_energy_on_column (const KT::MemberType& team, const int nlevs, const uview_1d& pseudo_density, @@ -348,7 +534,7 @@ compute_total_energy_on_column (const KT::MemberType& team, } KOKKOS_INLINE_FUNCTION -Real MassAndEnergyColumnConservationCheck:: +Real MassAndEnergyConservationCheck:: compute_energy_boundary_flux_on_column (const Real vapor_flux, const Real water_flux, const Real ice_flux, @@ -362,4 +548,17 @@ compute_energy_boundary_flux_on_column (const Real vapor_flux, return vapor_flux*(LatVap+LatIce) - (water_flux-ice_flux)*RHO_H2O*LatIce + heat_flux; } +Real MassAndEnergyConservationCheck::get_echeck() const{ + return echeck; +} + +Real MassAndEnergyConservationCheck::get_total_energy_before() const{ + return total_energy_before; +} + +Real MassAndEnergyConservationCheck::get_pb_fixer() const{ + return pb_fixer; +} + + } // namespace scream diff --git a/components/eamxx/src/share/property_checks/mass_and_energy_column_conservation_check.hpp b/components/eamxx/src/share/property_checks/mass_and_energy_conservation_check.hpp similarity index 86% rename from components/eamxx/src/share/property_checks/mass_and_energy_column_conservation_check.hpp rename to components/eamxx/src/share/property_checks/mass_and_energy_conservation_check.hpp index cd15cb4e0aa4..2dbae4d8e465 100644 --- a/components/eamxx/src/share/property_checks/mass_and_energy_column_conservation_check.hpp +++ b/components/eamxx/src/share/property_checks/mass_and_energy_conservation_check.hpp @@ -4,14 +4,16 @@ #include "share/property_checks/property_check.hpp" #include "share/grid/abstract_grid.hpp" #include "share/field/field.hpp" +#include "share/field/field_utils.hpp" -#include "ekat/kokkos/ekat_kokkos_utils.hpp" +#include +#include "ekat_comm.hpp" namespace scream { // This property check ensures that energy has been conserved. // It is a column-local check meant only for column independant processes. -class MassAndEnergyColumnConservationCheck: public PropertyCheck { +class MassAndEnergyConservationCheck: public PropertyCheck { using KT = KokkosTypes; using ExeSpaceUtils = ekat::ExeSpaceUtils; @@ -26,10 +28,13 @@ class MassAndEnergyColumnConservationCheck: public PropertyCheck { template using uview_2d = typename ekat::template Unmanaged >; + Field field_version_s1; + public: // Constructor - MassAndEnergyColumnConservationCheck (const std::shared_ptr& grid, + MassAndEnergyConservationCheck (const ekat::Comm& comm, + const std::shared_ptr& grid, const Real mass_error_tolerance, const Real energy_error_tolerance, const Field& pseudo_density_ptr, @@ -73,6 +78,12 @@ class MassAndEnergyColumnConservationCheck: public PropertyCheck { // in m_fields. void compute_current_energy (); + void global_fixer(const bool &); + + Real get_echeck() const; + Real get_total_energy_before() const; + Real get_pb_fixer() const; + // CUDA requires the parent fcn of a KOKKOS_LAMBDA to have public access #ifndef KOKKOS_ENABLE_CUDA protected: @@ -87,6 +98,11 @@ class MassAndEnergyColumnConservationCheck: public PropertyCheck { const uview_1d& qi, const uview_1d& qr); + KOKKOS_INLINE_FUNCTION + static Real compute_gas_mass_on_column (const KT::MemberType& team, + const int nlevs, + const uview_1d& pseudo_density); + KOKKOS_INLINE_FUNCTION static Real compute_mass_boundary_flux_on_column (const Real vapor_flux, const Real water_flux); @@ -109,9 +125,11 @@ class MassAndEnergyColumnConservationCheck: public PropertyCheck { const Real ice_flux, const Real heat_flux); + protected: std::shared_ptr m_grid; + ekat::Comm m_comm; std::map m_fields; int m_num_cols; @@ -120,10 +138,15 @@ class MassAndEnergyColumnConservationCheck: public PropertyCheck { Real m_mass_tol; Real m_energy_tol; + Real pb_fixer, echeck; + Real total_gas_mass_after, total_energy_before; + // Current value for total energy. These values // should be updated before a process is run. view_1d m_current_energy; view_1d m_current_mass; + + view_1d m_energy_change; }; // class EnergyConservationCheck } // namespace scream diff --git a/components/eamxx/src/share/util/eamxx_repro_sum_mod.F90 b/components/eamxx/src/share/util/eamxx_repro_sum_mod.F90 new file mode 100644 index 000000000000..48bd9e6b2ac3 --- /dev/null +++ b/components/eamxx/src/share/util/eamxx_repro_sum_mod.F90 @@ -0,0 +1,22 @@ +module eamxx_repro_sum_mod + + implicit none + +contains + + subroutine eamxx_repro_sum(send, recv, nlocal, nfld, comm) bind(c) + use iso_c_binding, only: c_int, c_double + + use shr_reprosum_mod, only: repro_sum => shr_reprosum_calc + + integer(kind=c_int), value, intent(in) :: nlocal, nfld, comm + real(kind=c_double), intent(in) :: send(nlocal,nfld) + real(kind=c_double), intent(out) :: recv(nfld) + + call repro_sum(send, recv, nlocal, nlocal, nfld, commid=comm) + end subroutine eamxx_repro_sum + +end module eamxx_repro_sum_mod + + + diff --git a/components/eamxx/src/share/util/eamxx_utils.hpp b/components/eamxx/src/share/util/eamxx_utils.hpp index 2b30a1c471ac..a8776b17c324 100644 --- a/components/eamxx/src/share/util/eamxx_utils.hpp +++ b/components/eamxx/src/share/util/eamxx_utils.hpp @@ -443,6 +443,10 @@ struct DefaultMetadata { } }; +extern "C" +void eamxx_repro_sum(const Real* send, Real* recv, + Int nlocal, Int nfld, Int fcomm); + } // namespace scream #endif // SCREAM_UTILS_HPP diff --git a/share/util/shr_reprosum_mod.F90 b/share/util/shr_reprosum_mod.F90 index a0bd5905a2f8..e0a3bf954a6b 100644 --- a/share/util/shr_reprosum_mod.F90 +++ b/share/util/shr_reprosum_mod.F90 @@ -45,7 +45,9 @@ module shr_reprosum_mod shr_infnan_nan, & shr_infnan_isnan, shr_infnan_isinf, & shr_infnan_isposinf, shr_infnan_isneginf +#ifndef EAMXX_STANDALONE use perf_mod +#endif ! Import MPI fcns/types use mpi @@ -92,6 +94,23 @@ module shr_reprosum_mod logical :: repro_sum_allow_infnan = .false. +#ifdef EAMXX_STANDALONE + ! Declare the C function interface + interface + subroutine shr_reprosumx86_fix_start(arg) bind(c) + use iso_c_binding + integer, intent(out) :: arg + end subroutine shr_reprosumx86_fix_start + end interface + + interface + subroutine shr_reprosumx86_fix_end(arg) bind(c) + use iso_c_binding + integer, intent(in) :: arg + end subroutine shr_reprosumx86_fix_end + end interface +#endif + CONTAINS ! @@ -604,7 +623,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & else mpi_comm = MPI_COMM_WORLD endif +#ifndef EAMXX_STANDALONE call t_barrierf('sync_repro_sum',mpi_comm) +#endif ! Check whether should abort if input contains NaNs or INFs abort_inf_nan = .not. repro_sum_allow_infnan @@ -617,7 +638,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & abort_inf_nan = .true. #endif +#ifndef EAMXX_STANDALONE call t_startf('shr_reprosum_INF_NaN_Chk') +#endif ! Initialize flags to indicate that no NaNs or INFs are present in the input data inf_nan_gchecks = .false. @@ -656,12 +679,15 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & inf_nan_lchecks(2,ifld) = any(shr_infnan_isposinf(arr(:,ifld))) inf_nan_lchecks(3,ifld) = any(shr_infnan_isneginf(arr(:,ifld))) end do - +#ifndef EAMXX_STANDALONE call t_startf("repro_sum_allr_lor") +#endif call mpi_allreduce (inf_nan_lchecks, inf_nan_gchecks, 3*nflds, & MPI_LOGICAL, MPI_LOR, mpi_comm, ierr) gbl_lor_red = 1 +#ifndef EAMXX_STANDALONE call t_stopf("repro_sum_allr_lor") +#endif do ifld=1,nflds arr_gsum_infnan(ifld) = any(inf_nan_gchecks(:,ifld)) @@ -670,7 +696,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & endif +#ifndef EAMXX_STANDALONE call t_stopf('shr_reprosum_INF_NaN_Chk') +#endif ! Check whether should use shr_reprosum_ddpdd algorithm use_ddpdd_sum = repro_sum_use_ddpdd @@ -685,17 +713,23 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & if ( use_ddpdd_sum ) then +#ifndef EAMXX_STANDALONE call t_startf('shr_reprosum_ddpdd') +#endif call shr_reprosum_ddpdd(arr, arr_gsum, nsummands, dsummands, & nflds, mpi_comm) repro_sum_fast = 1 +#ifndef EAMXX_STANDALONE call t_stopf('shr_reprosum_ddpdd') +#endif else +#ifndef EAMXX_STANDALONE call t_startf('shr_reprosum_int') +#endif ! Get number of MPI tasks call mpi_comm_size(mpi_comm, tasks, ierr) @@ -733,7 +767,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & ! Determine maximum number of summands in local phases of the ! algorithm +#ifndef EAMXX_STANDALONE call t_startf("repro_sum_allr_max") +#endif if ( present(gbl_max_nsummands) ) then if (gbl_max_nsummands < 1) then call mpi_allreduce (nsummands, max_nsummands, 1, & @@ -747,7 +783,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & MPI_INTEGER, MPI_MAX, mpi_comm, ierr) gbl_max_red = 1 endif +#ifndef EAMXX_STANDALONE call t_stopf("repro_sum_allr_max") +#endif ! Determine maximum shift. Shift needs to be small enough that summation, ! in absolute value, does not exceed maximum value representable by i8. @@ -851,7 +889,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & !$omp default(shared) & !$omp private(ithread, ifld, isum, arr_exp, arr_exp_tlmin, arr_exp_tlmax) do ithread=1,omp_nthreads +#ifndef EAMXX_STANDALONE call t_startf('repro_sum_loopa') +#endif do ifld=1,nflds arr_exp_tlmin = MAXEXPONENT(1.0_r8) arr_exp_tlmax = MINEXPONENT(1.0_r8) @@ -867,7 +907,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & arr_tlmin_exp(ifld,ithread) = arr_exp_tlmin arr_tlmax_exp(ifld,ithread) = arr_exp_tlmax end do +#ifndef EAMXX_STANDALONE call t_stopf('repro_sum_loopa') +#endif end do do ifld=1,nflds @@ -879,10 +921,14 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & arr_lextremes(0,:) = -nsummands arr_lextremes(1:nflds,1) = -arr_lmax_exp(:) arr_lextremes(1:nflds,2) = arr_lmin_exp(:) +#ifndef EAMXX_STANDALONE call t_startf("repro_sum_allr_minmax") +#endif call mpi_allreduce (arr_lextremes, arr_gextremes, 2*(nflds+1), & MPI_INTEGER, MPI_MIN, mpi_comm, ierr) +#ifndef EAMXX_STANDALONE call t_stopf("repro_sum_allr_minmax") +#endif max_nsummands = -arr_gextremes(0,1) arr_gmax_exp(:) = -arr_gextremes(1:nflds,1) arr_gmin_exp(:) = arr_gextremes(1:nflds,2) @@ -999,16 +1045,21 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & endif +#ifndef EAMXX_STANDALONE call t_stopf('shr_reprosum_int') +#endif endif ! Compare integer vector and floating point results if ( present(rel_diff) ) then if (shr_reprosum_reldiffmax >= 0.0_r8) then - +#ifndef EAMXX_STANDALONE call t_barrierf('sync_nonrepro_sum',mpi_comm) +#endif +#ifndef EAMXX_STANDALONE call t_startf('nonrepro_sum') +#endif ! Record statistic nonrepro_sum = 1 ! Compute nonreproducible sum @@ -1024,12 +1075,18 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & endif end do +#ifndef EAMXX_STANDALONE call t_startf("nonrepro_sum_allr_r8") +#endif call mpi_allreduce (arr_lsum, arr_gsum_fast, nflds, & MPI_REAL8, MPI_SUM, mpi_comm, ierr) +#ifndef EAMXX_STANDALONE call t_stopf("nonrepro_sum_allr_r8") +#endif +#ifndef EAMXX_STANDALONE call t_stopf('nonrepro_sum') +#endif ! Determine differences !$omp parallel do & @@ -1271,7 +1328,9 @@ subroutine shr_reprosum_int (arr, arr_gsum, nsummands, dsummands, nflds, & !$omp private(ithread, ifld, ioffset, isum, arr_frac, arr_exp, & !$omp arr_shift, ilevel, i8_arr_level, arr_remainder, RX_8, IX_8) do ithread=1,omp_nthreads +#ifndef EAMXX_STANDALONE call t_startf('repro_sum_loopb') +#endif do ifld=1,nflds ioffset = offset(ifld) @@ -1414,7 +1473,9 @@ subroutine shr_reprosum_int (arr, arr_gsum, nsummands, dsummands, nflds, & endif enddo enddo +#ifndef EAMXX_STANDALONE call t_stopf('repro_sum_loopb') +#endif enddo ! Sum contributions from different threads @@ -1442,18 +1503,28 @@ subroutine shr_reprosum_int (arr, arr_gsum, nsummands, dsummands, nflds, & ! Sum integer vector element-wise #if ( defined noI8 ) ! Workaround for when shr_kind_i8 is not supported. +#ifndef EAMXX_STANDALONE call t_startf("repro_sum_allr_i4") +#endif call mpi_allreduce (i8_arr_lsum_level, i8_arr_gsum_level, & veclth, MPI_INTEGER, MPI_SUM, mpi_comm, ierr) +#ifndef EAMXX_STANDALONE call t_stopf("repro_sum_allr_i4") +#endif #else +#ifndef EAMXX_STANDALONE call t_startf("repro_sum_allr_i8") +#endif call mpi_allreduce (i8_arr_lsum_level, i8_arr_gsum_level, & veclth, MPI_INTEGER8, MPI_SUM, mpi_comm, ierr) +#ifndef EAMXX_STANDALONE call t_stopf("repro_sum_allr_i8") #endif +#endif +#ifndef EAMXX_STANDALONE call t_startf('repro_sum_finalsum') +#endif ! Construct global sum from integer vector representation: ! 1) arr_max_shift is the shift applied to fraction(arr_gmax) . ! When shifting back, need to 'add back in' the true arr_gmax exponent. @@ -1827,7 +1898,9 @@ subroutine shr_reprosum_int (arr, arr_gsum, nsummands, dsummands, nflds, & endif enddo +#ifndef EAMXX_STANDALONE call t_stopf('repro_sum_finalsum') +#endif end subroutine shr_reprosum_int @@ -1990,10 +2063,14 @@ subroutine shr_reprosum_ddpdd (arr, arr_gsum, nsummands, dsummands, & enddo +#ifndef EAMXX_STANDALONE call t_startf("repro_sum_allr_c16") +#endif call mpi_allreduce (arr_lsum_dd, arr_gsum_dd, nflds, & MPI_COMPLEX16, mpi_sumdd, mpi_comm, ierr) +#ifndef EAMXX_STANDALONE call t_stopf("repro_sum_allr_c16") +#endif do ifld=1,nflds arr_gsum(ifld) = real(arr_gsum_dd(ifld))