EAMxx: dycore energy fixer

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",
             )
         },

components/eamxx/cime_config/namelist_defaults_eamxx.xml

@@ -193,6 +193,10 @@ be lost if SCREAM_HACK_XML is not enabled.
       <moisture>moist</moisture>
       <!-- Frequency in physics steps to output a global hash over the dycore's
            in-fields. <= 0 disables hashing. -->
+      <enable_column_conservation_checks>false</enable_column_conservation_checks>
+      <enable_energy_fixer type="logical" doc="Turn on energy fixer for dycore (fixes dycore, pressure adjustment, PD coupling)">true</enable_energy_fixer>
+      <!-- true means 2 more calls to AllReduce after dycore -->
+      <enable_energy_fixer_debug_info type="logical" doc="Prints debug info that energy fixer indeed fixed energy (cost of an MPI reduce)">false</enable_energy_fixer_debug_info>
       <bfb_hash type="integer">18</bfb_hash>
     </homme>
 

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"

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)
       $<$<AND:$<COMPILE_LANGUAGE:Fortran>,$<CXX_COMPILER_ID:GNU>>:CPRGNU>
       $<$<AND:$<COMPILE_LANGUAGE:Fortran>,$<CXX_COMPILER_ID:Intel>>:CPRINTEL>
       $<$<AND:$<COMPILE_LANGUAGE:Fortran>,$<CXX_COMPILER_ID:Clang>>:CPRCRAY>
-      $<$<AND:$<COMPILE_LANGUAGE:Fortran>,$<CXX_COMPILER_ID:clang>>:CPRCRAY>)
+      $<$<AND:$<COMPILE_LANGUAGE:Fortran>,$<CXX_COMPILER_ID:clang>>:CPRCRAY>
+      EAMXX_STANDALONE)
 
 
     if (${CMAKE_SYSTEM} MATCHES "Linux")

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} <CASE_DIR> [<param>=<val>] ...
+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__)

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<MassAndEnergyColumnConservationCheck>(phys_grid,
+    std::make_shared<MassAndEnergyConservationCheck>(m_atm_comm,phys_grid,
                                                            mass_error_tol, energy_error_tol,
                                                            pseudo_density, ps, phis,
                                                            horiz_winds, T_mid, qv,

components/eamxx/src/dynamics/homme/eamxx_homme_process_interface.cpp

@@ -254,7 +254,18 @@ void HommeDynamics::set_grids (const std::shared_ptr<const GridsManager> 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<Computed>("vapor_flux", pg_scalar2d, kg/(m2*s), grid_name);
+    add_field<Computed>("water_flux", pg_scalar2d, m/s,     grid_name);
+    add_field<Computed>("ice_flux",   pg_scalar2d, m/s,     grid_name);
+    add_field<Computed>("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::

components/eamxx/src/share/CMakeLists.txt

@@ -34,14 +34,15 @@ 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
   util/eamxx_time_stamp.cpp
   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()