Merge pull request #260 from eagles-project/oscar/setext

Oscar/setext

src/mam4xx/CMakeLists.txt

@@ -43,6 +43,7 @@ install(FILES
         mo_photo.hpp
         lin_strat_chem.hpp
         mo_chm_diags.hpp
+        mo_setext.hpp
         DESTINATION include/mam4xx)
 
 add_library(mam4xx aero_modes.cpp)

src/mam4xx/mam4.hpp

@@ -24,6 +24,7 @@
 #include <mam4xx/mam4_types.hpp>
 #include <mam4xx/mo_chm_diags.hpp>
 #include <mam4xx/mo_photo.hpp>
+#include <mam4xx/mo_setext.hpp>
 #include <mam4xx/ndrop.hpp>
 #include <mam4xx/nucleate_ice.hpp>
 #include <mam4xx/nucleation.hpp>

src/mam4xx/mo_setext.hpp

@@ -0,0 +1,142 @@
+#ifndef MAM4XX_MO_SETEXT_HPP
+#define MAM4XX_MO_SETEXT_HPP
+
+#include <haero/math.hpp>
+#include <mam4xx/aero_config.hpp>
+#include <mam4xx/mam4_types.hpp>
+#include <mam4xx/utils.hpp>
+
+namespace mam4 {
+
+namespace mo_setext {
+using View1D = DeviceType::view_1d<Real>;
+using View2D = DeviceType::view_2d<Real>;
+using View3D = DeviceType::view_3d<Real>;
+constexpr int pver = mam4::nlev;
+constexpr int extfrc_cnt = 9;
+constexpr int extcnt = 9; //, & ! number of species with external forcing
+
+struct Forcing {
+  int frc_ndx;
+  bool file_alt_data;
+  View2D fields_data;
+  int nsectors;
+};
+
+KOKKOS_INLINE_FUNCTION
+void extfrc_set(const Forcing *forcings, const View2D &frcing) {
+
+  /*--------------------------------------------------------
+   ... form the external forcing
+  --------------------------------------------------------*/
+  // param[in] forcings(extcnt) array with a list of Forcing object.
+  // @param[out] frcing(ncol,pver,extcnt)   insitu forcings [molec/cm^3/s]
+  // Note: we do not need zint to compute frcing
+  // const ColumnView &zint
+
+  constexpr Real zero = 0.0;
+
+  // frcing(:,:,:) = zero;
+  if (extfrc_cnt < 1 || extcnt < 1) {
+    return;
+  }
+
+  /*--------------------------------------------------------
+  ! ... set non-zero forcings
+  !--------------------------------------------------------*/
+
+  for (int mm = 0; mm < extfrc_cnt; ++mm) {
+    // Fortran to C++ indexing
+    auto forcing_mm = forcings[mm];
+    const int nn = forcing_mm.frc_ndx - 1;
+    for (int kk = 0; kk < pver; ++kk) {
+      frcing(kk, nn) = zero;
+    } // k
+
+    for (int isec = 0; isec < forcing_mm.nsectors; ++isec) {
+      if (forcing_mm.file_alt_data) {
+        for (int kk = 0; kk < pver; ++kk) {
+          // frcing(:ncol,:,nn) = frcing(:ncol,:,nn) + &
+          // forcings(mm)%fields(isec)%data(:ncol,pver:1:-1,lchnk)
+          frcing(kk, nn) += forcing_mm.fields_data(isec, pver - 1 - kk);
+        } // kk
+      } else {
+        // // forcings(mm)%fields(isec)%data(:ncol,:,lchnk)
+        for (int kk = 0; kk < pver; ++kk) {
+          frcing(kk, nn) += forcing_mm.fields_data(isec, kk);
+        }
+      }
+    } // isec
+
+    // I did not include variables that are involved in the outfld.
+    //
+    //  xfcname = trim(forcings(mm)%species)//'_XFRC'
+    // call outfld( xfcname, frcing(:ncol,:,nn), ncol, lchnk )
+
+    // frcing_col(:ncol) = 0._r8
+    // do kk = 1,pver
+    //    frcing_col(:ncol) = frcing_col(:ncol) + &
+    //                  frcing(:ncol,kk,nn)*(zint(:ncol,kk)-zint(:ncol,kk+1))*km_to_cm
+    // enddo
+    // xfcname = trim(forcings(mm)%species)//'_CLXF'
+    // call outfld( xfcname, frcing_col(:ncol), ncol, lchnk )
+
+  } // end mm
+
+} // extfrc_set
+
+KOKKOS_INLINE_FUNCTION
+void setext(const Forcing *forcings,
+            const View2D &extfrc) // ! out
+{
+  // Note: we do not need setext.
+  // Because I removed diagnostic variables(i.e., output variable in outfld),
+  // setext and extfrc_set are equivalent.
+  /*--------------------------------------------------------
+       ... for this latitude slice:
+           - form the production from datasets
+           - form the nox (xnox) production from lighing
+           - form the nox (xnox) production from airplanes
+           - form the co production from airplanes
+  --------------------------------------------------------*/
+
+  // param[in] forcings(extcnt) array with a list of Forcing object.
+  // @param[out]  extfrc(ncol,pver,extcnt)    ! the "extraneous" forcing
+
+  /*--------------------------------------------------------
+    !     ... local variables
+    !--------------------------------------------------------
+    ! variables for output. in current MAM4 they are not calculated and are
+    assigned zero real(r8), dimension(ncol,pver) :: no_lgt, no_air, co_air */
+
+  /*--------------------------------------------------------
+  !     ... set frcing from datasets
+  !--------------------------------------------------------*/
+  extfrc_set(forcings,
+             extfrc); // out
+
+  /*--------------------------------------------------------
+       ... set nox production from lighting
+           note: from ground to cloud top production is c shaped
+
+   FORTRAN refactor: nox is not included in current MAM4
+   the related code are removed but outfld is kept for BFB testing
+  --------------------------------------------------------*/
+  // no_lgt(:,:) = 0._r8
+  // call outfld( 'NO_Lightning', no_lgt(:ncol,:), ncol, lchnk )
+
+  // ! FORTRAN refactor: in the subroutine airpl_set, has_airpl_src is false,
+  // ! the subroutine only has two outfld calls that output zero
+  // ! remove the subroutine call and move out the zero outfld
+  // no_air(:,:) = 0._r8
+  // co_air(:,:) = 0._r8
+  // call outfld( 'NO_Aircraft',  no_air(:ncol,:), ncol, lchnk )
+  // call outfld( 'CO_Aircraft',  co_air(:ncol,:), ncol, lchnk )
+
+} // setext
+
+} // namespace mo_setext
+
+} // end namespace mam4
+
+#endif

src/validation/CMakeLists.txt

@@ -31,3 +31,4 @@ add_subdirectory(water_uptake)
 add_subdirectory(mo_photo)
 add_subdirectory(lin_strat_chem)
 add_subdirectory(mo_chm_diags)
+add_subdirectory(mo_setext)

src/validation/mo_setext/CMakeLists.txt

@@ -0,0 +1,55 @@
+set(SETEXT_VALIDATION_DIR ${MAM_X_VALIDATION_DIR}/mo_setext)
+set(SETEXT_VALIDATION_SCRIPTS_DIR ${MAM_X_VALIDATION_DIR}/scripts)
+
+
+# These subdirectories contain Skywalker drivers for MAM4 parameterizations.
+# Include directory for .mod files.
+
+include_directories(${PROJECT_BINARY_DIR}/validation)
+
+# We use a single driver for all setext-related parameterizations.
+add_executable(setext_driver
+               setext_driver.cpp
+               extfrc_set.cpp
+               )
+
+target_link_libraries(setext_driver skywalker;validation;${HAERO_LIBRARIES})
+
+# Copy some Python scripts from mam_x_validation to our binary directory.
+foreach(script
+        compare_mam4xx_mam4.py)
+  configure_file(
+    ${SETEXT_VALIDATION_SCRIPTS_DIR}/${script}
+    ${CMAKE_CURRENT_BINARY_DIR}/${script}
+    COPYONLY
+  )
+endforeach()
+# stand_calc_sum_wght_ts_355
+# Run the driver in several configurations to produce datasets.
+set(TEST_LIST
+    extfrc_set_ts_355
+    )
+# # matching the tests and errors, just for convenience
+
+set(DEFAULT_TOL 1e-13)
+set(ERROR_THRESHOLDS
+9e-10
+   )
+foreach(input tol IN ZIP_LISTS TEST_LIST ERROR_THRESHOLDS)
+  # copy the baseline file into place.
+
+  configure_file(
+    ${SETEXT_VALIDATION_DIR}/mam_${input}.py
+    ${CMAKE_CURRENT_BINARY_DIR}/mam_${input}.py
+    COPYONLY
+  )
+
+  # add a test to run the skywalker driver
+  add_test(run_${input} setext_driver ${SETEXT_VALIDATION_DIR}/${input}.yaml)
+
+  # add a test to validate mam4xx's results against the baseline.
+  # Select a threshold error slightly bigger than the largest relative error for the threshold error.
+  # compare_mam4xx_mam4.py <module1.py> <module2.py> <check_norms> <threshold error>
+  add_test(validate_${input} python3 compare_mam4xx_mam4.py mam4xx_${input}.py mam_${input}.py True ${tol})
+  set_tests_properties(validate_${input} PROPERTIES DEPENDS run_${input})
+endforeach()

src/validation/mo_setext/extfrc_set.cpp

@@ -0,0 +1,64 @@
+// mam4xx: Copyright (c) 2022,
+// Battelle Memorial Institute and
+// National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <mam4xx/mam4.hpp>
+
+#include <mam4xx/aero_config.hpp>
+#include <skywalker.hpp>
+#include <validation.hpp>
+
+using namespace skywalker;
+using namespace mam4;
+using namespace haero;
+using namespace mo_setext;
+void extfrc_set(Ensemble *ensemble) {
+  ensemble->process([=](const Input &input, Output &output) {
+    using View1DHost = typename HostType::view_1d<Real>;
+
+    Forcing forcings[extfrc_cnt];
+    for (int i = 1; i <= extfrc_cnt; ++i) {
+
+      Forcing forcing_mm;
+
+      forcing_mm.frc_ndx =
+          int(input.get_array("forcings" + std::to_string(i) + "_frc_ndx")[0]);
+      forcing_mm.nsectors =
+          int(input.get_array("forcings" + std::to_string(i) + "_nsectors")[0]);
+      forcing_mm.file_alt_data = int(input.get_array(
+          "forcings" + std::to_string(i) + "_file_alt_data")[0]);
+
+      forcing_mm.fields_data = View2D("data", forcing_mm.nsectors, pver);
+
+      for (int isec = 1; isec <= forcing_mm.nsectors; ++isec) {
+
+        auto label = "forcings" + std::to_string(i) + "_fields" +
+                     std::to_string(isec) + "_data";
+        const auto data1 = input.get_array(label);
+
+        View1DHost forcings_fields_data_host =
+            View1DHost((Real *)data1.data(), data1.size());
+
+        const View1D data_isec =
+            Kokkos::subview(forcing_mm.fields_data, isec - 1, Kokkos::ALL());
+        Kokkos::deep_copy(data_isec, forcings_fields_data_host);
+      } // isec
+
+      forcings[i - 1] = forcing_mm;
+    }
+
+    View2D frcing("frcing", pver, extcnt);
+    const int ncol = 1;
+    auto team_policy = ThreadTeamPolicy(ncol, 1u);
+
+    Kokkos::parallel_for(
+        team_policy, KOKKOS_LAMBDA(const ThreadTeam &team) {
+          extfrc_set(forcings, frcing);
+        });
+    std::vector<Real> frcing_out(pver * extcnt, 0.0);
+    mam4::validation::convert_2d_view_device_to_1d_vector(frcing, frcing_out);
+
+    output.set("frcing", frcing_out);
+  });
+}

src/validation/mo_setext/setext_driver.cpp

@@ -0,0 +1,68 @@
+// mam4xx: Copyright (c) 2022,
+// Battelle Memorial Institute and
+// National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <mam4xx/mo_setext.hpp>
+
+#include <iostream>
+#include <skywalker.hpp>
+#include <validation.hpp>
+
+void usage() {
+  std::cerr << "setext_driver: a Skywalker driver for validating the "
+               "MAM4 rename parameterizations."
+            << std::endl;
+  std::cerr << "setext_driver: usage:" << std::endl;
+  std::cerr << "setext_driver <input.yaml>" << std::endl;
+  exit(0);
+}
+
+using namespace skywalker;
+using namespace mam4;
+
+// Parameterizations used by the setext() process.
+void extfrc_set(Ensemble *ensemble);
+
+int main(int argc, char **argv) {
+  if (argc == 1) {
+    usage();
+  }
+  validation::initialize(argc, argv);
+  std::string input_file = argv[1];
+  std::string output_file = validation::output_name(input_file);
+  std::cout << argv[0] << ": reading " << input_file << std::endl;
+
+  // Load the ensemble. Any error encountered is fatal.
+  Ensemble *ensemble = skywalker::load_ensemble(input_file, "mam4xx");
+
+  // the settings.
+  Settings settings = ensemble->settings();
+  if (!settings.has("function")) {
+    std::cerr << "No function specified in mam4xx.settings!" << std::endl;
+    exit(1);
+  }
+
+  // Dispatch to the requested function.
+  auto func_name = settings.get("function");
+  try {
+    if (func_name == "extfrc_set") {
+      extfrc_set(ensemble);
+    } else {
+      std::cerr << "Error: Function name '" << func_name
+                << "' does not have an implemented test!" << std::endl;
+      exit(1);
+    }
+
+  } catch (std::exception &e) {
+    std::cerr << argv[0] << ": Error: " << e.what() << std::endl;
+  }
+
+  // Write out a Python module.
+  std::cout << argv[0] << ": writing " << output_file << std::endl;
+  ensemble->write(output_file);
+
+  // Clean up.
+  delete ensemble;
+  validation::finalize();
+}