diff --git a/CMakeLists.txt b/CMakeLists.txt index 3bc54a375..75a29f4a5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,7 @@ project(nextsim_dg) option(WITH_THREADS "Build with support for openmp" OFF) option(ENABLE_MPI "Enable distributed-memory parallelization with MPI" OFF) option(ENABLE_XIOS "Enable XIOS library for IO" OFF) +option(ENABLE_OASIS "Enable the OASIS interface" OFF) option(BUILD_TESTS "Build the tests" ON) set(DynamicsType "DG2" @@ -145,6 +146,23 @@ if(ENABLE_XIOS) ) endif() +# Do we want to include the OASIS interface? It'll only work if we have MPI +if (ENABLE_OASIS) + if (ENABLE_MPI) + message(STATUS "Building with OASIS support") + find_package(OASIS REQUIRED) + pkg_search_module(NETCDF REQUIRED netcdf-fortran) + target_compile_definitions(nextsimlib PUBLIC USE_OASIS) + + target_link_directories(nextsimlib PUBLIC ${NETCDF_LIBRARY_DIRS}) + target_link_directories(nextsimlib PUBLIC ${OASIS_LIBRARIES}) + target_link_libraries(nextsimlib PUBLIC oasis.cbind mct mpeu psmile.MPI1 scrip netcdff netcdf) + target_include_directories(nextsimlib PUBLIC ${OASIS_INCLUDES}) + else () + message(FATAL_ERROR "Cannot build with OASIS support, because MPI is not enabled" .) + endif () +endif () + # OpenMP if(WITH_THREADS) find_package(OpenMP REQUIRED) diff --git a/cmake/FindOASIS.cmake b/cmake/FindOASIS.cmake new file mode 100644 index 000000000..ecf2d3478 --- /dev/null +++ b/cmake/FindOASIS.cmake @@ -0,0 +1,21 @@ +# Find oasis +# +# Please pass the -DOASIS_DIR variable to cmake (location of the OASIS libraries). + +find_library (OASIS_LIBRARIES NAMES mct mpeu oasis psmile scrip HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH) + +get_filename_component (OASIS_LIBRARIES "${OASIS_LIBRARIES}" PATH) + +set (OASIS_DIR "${OASIS_LIBRARIES}/../") +cmake_path(NORMAL_PATH OASIS_DIR) + +find_path (OASIS_INCLUDES NAMES oasis_c.h HINTS ${OASIS_DIR}/include) + +# handle the QUIETLY and REQUIRED arguments and set OASIS_FOUND to TRUE if all listed variables are TRUE +include (FindPackageHandleStandardArgs) +find_package_handle_standard_args(OASIS DEFAULT_MSG + OASIS_LIBRARIES OASIS_INCLUDES) + +# message (STATUS "OASIS_LIBRARIES: ${OASIS_LIBRARIES}") +# message (STATUS "OASIS_INCLUDES: ${OASIS_INCLUDES}") +# message (STATUS "OASIS_DIR: ${OASIS_DIR}") diff --git a/core/src/Model.cpp b/core/src/Model.cpp index 63bf2f9e6..17cd46963 100644 --- a/core/src/Model.cpp +++ b/core/src/Model.cpp @@ -1,6 +1,6 @@ /*! * @file Model.cpp - * @date 12 Aug 2021 + * @date 09 Sep 2024 * @author Tim Spain * @author Kacper Kornet */ @@ -47,6 +47,9 @@ static std::map modelConfigKeyMap = { { Model::MISSINGVALUE_KEY, "model.missing_value" }, { Model::RESTARTPERIOD_KEY, "model.restart_period" }, { Model::RESTARTOUTFILE_KEY, "model.restart_file" }, +#ifdef USE_OASIS + { Model::WRITEOASISGRID_KEY, "oasis.write_grid" }, +#endif }; // Merge the configuration from ModelConfig into the Model keyMap. @@ -136,10 +139,16 @@ void Model::configure() // Get the coordinates from the ModelState for persistence m_etadata.extractCoordinates(initialState); +#ifdef USE_OASIS + const bool writeOasisGrid = Configured::getConfiguration(keyMap.at(WRITEOASISGRID_KEY), false); + m_etadata.initOasis(writeOasisGrid); +#endif + modelStep.setData(pData); modelStep.setMetadata(m_etadata); modelStep.setRestartDetails(restartPeriod, finalFileName); pData.setData(initialState.data); + pData.setMetadata(m_etadata); } ConfigMap Model::getConfig() const diff --git a/core/src/ModelMetadata.cpp b/core/src/ModelMetadata.cpp index d30597ec6..39340e08f 100644 --- a/core/src/ModelMetadata.cpp +++ b/core/src/ModelMetadata.cpp @@ -1,7 +1,7 @@ /*! * @file ModelMetadata.cpp * - * @date 21 August 2024 + * @date 10 Sep 2024 * @author Tim Spain */ @@ -16,6 +16,9 @@ #include #include #endif +#ifdef USE_OASIS +#include +#endif namespace Nextsim { @@ -102,4 +105,70 @@ ModelState& ModelMetadata::affixCoordinates(ModelState& state) const } return state; } + +#ifdef USE_OASIS +void ModelMetadata::initOasis(const bool writeOasisGrid) +{ + // Set the partitioning + /* From the manual: "[ig_paral is a] vector of integers describing the local grid partition + * in the global index space; has a different expression depending on the type of the + * partition; in OASIS3-MCT, 5 types of partition are supported: Serial (no partition), + * Apple, Box, Orange, and Points" - it looks like we should use "Box", so partInfo[0] = 2 + * (aka. ig_paral). + * + * #Box partition# + * Each partition is a rectangular region of the global domain, described by the global + * offset of its upper left corner, and its local extents in the X and Y dimensions. The + * global extent in the X dimension must also be given. In this case, we have ig_paral(1:5): + * - ig_paral(1) = 2 (indicates a Box partition) + * - ig_paral(2) = the upper left corner global offset + * - ig paral(3) = the local extent in x + * - ig_paral(4) = the local extent in y + * - ig_paral(5) = the global extent in x. + * + * metdatata contains: localCornerX, localCornerY, localExtentX, localExtentY, + * globalExtentX, globalExtentY; + */ + // TODO: The contents of metadata is not certain! + const int offset = globalExtentX * localCornerY + localCornerX; + const std::vector partInfo + = { OASIS_Box, offset, localExtentX, localExtentY, globalExtentX }; + + const int globalSize = globalExtentX * globalExtentY; + const std::string compName = "nextsim"; // Not useful for any setups we have in mind + OASIS_CHECK_ERR(oasis_c_def_partition( + &OASISPartitionId, OASIS_Box_Params, &partInfo[0], globalSize, compName.c_str())); + + // TODO: Writing out grid information should be possible, but optional + if (writeOasisGrid) { + /* This needs to be figured out, but it's not a priority. Grid writing is + * not necessary for the type of coupling we'll start with. + + const std::string gridName = "nxts"; + + int flag = 1; + OASIS_CHECK_ERR(oasis_c_start_grids_writing(&flag)); + + OASIS_CHECK_ERR(oasis_c_write_grid( + gridName.c_str(), nx, ny, nx_loc, ny_loc, lon, lat, OASISPartitionId)); + OASIS_CHECK_ERR(oasis_c_write_corner( + gridName.c_str(), nx, ny, nx_loc, ny_loc, clo, cla, OASISPartitionId)); + OASIS_CHECK_ERR(oasis_c_write_area( + gridName.c_str(), nx, ny, nx_loc, ny_loc, area, OASISPartitionId)); + OASIS_CHECK_ERR(oasis_c_write_mask( + gridName.c_str(), nx, ny, nx_loc, ny_loc, angle, OASISPartitionId)); + + std::string companion = "land area fraction"; + OASIS_CHECK_ERR(oasis_c_write_frac( + gridName.c_str(), nx, ny, nx_loc, ny_loc, mask, OASISPartitionId), + companion.c_str()); + companion = "land sea mask"; + OASIS_CHECK_ERR(oasis_c_write_mask( + gridName.c_str(), nx, ny, nx_loc, ny_loc, mask, OASISPartitionId), + companion.c_str()); + */ + } +} +#endif + } /* namespace Nextsim */ diff --git a/core/src/PrognosticData.cpp b/core/src/PrognosticData.cpp index 9ac0b146e..3faa8f390 100644 --- a/core/src/PrognosticData.cpp +++ b/core/src/PrognosticData.cpp @@ -1,7 +1,7 @@ /*! * @file PrognosticData.cpp * - * @date 1 Jul 2024 + * @date 09 Sep 2024 * @author Tim Spain * @author Einar Ólason */ @@ -12,6 +12,10 @@ #include "include/Module.hpp" #include "include/gridNames.hpp" +#ifdef USE_OASIS +#include +#endif + namespace Nextsim { PrognosticData::PrognosticData() @@ -74,6 +78,17 @@ void PrognosticData::setData(const ModelState::DataMap& ms) iceGrowth.setData(ms); } +void PrognosticData::setMetadata(const Nextsim::ModelMetadata& metadata) +{ + pAtmBdy->setMetadata(metadata); + pOcnBdy->setMetadata(metadata); + +#ifdef USE_OASIS + // OASIS finalising definition - can only be called once + OASIS_CHECK_ERR(oasis_c_enddef()); +#endif +} + void PrognosticData::update(const TimestepTime& tst) { ModelArrayRef ticeUpd(getStore()); diff --git a/core/src/include/Model.hpp b/core/src/include/Model.hpp index 83eaab5e6..b71050f10 100644 --- a/core/src/include/Model.hpp +++ b/core/src/include/Model.hpp @@ -1,6 +1,6 @@ /*! * @file Model.hpp - * @date 12 Aug 2021 + * @date 09 Sep 2024 * @author Tim Spain * @author Kacper Kornet */ @@ -49,6 +49,9 @@ class Model : public Configured { // Other Model configuration keys, not to be written to the restart file. RESTARTPERIOD_KEY, RESTARTOUTFILE_KEY, +#ifdef USE_OASIS + WRITEOASISGRID_KEY, +#endif }; ConfigMap getConfig() const; diff --git a/core/src/include/ModelMetadata.hpp b/core/src/include/ModelMetadata.hpp index ad4956c7a..043d941a8 100644 --- a/core/src/include/ModelMetadata.hpp +++ b/core/src/include/ModelMetadata.hpp @@ -1,7 +1,7 @@ /*! * @file ModelMetadata.hpp * - * @date Jun 29, 2022 + * @date 10 Sep 2024 * @author Tim Spain */ @@ -42,6 +42,10 @@ class ModelMetadata { ModelMetadata() = default; #endif +#ifdef USE_OASIS + int OASISPartitionId; +#endif + /*! * @brief Sets the initial or current model time * @@ -99,6 +103,10 @@ class ModelMetadata { int localCornerX, localCornerY, localExtentX, localExtentY, globalExtentX, globalExtentY; #endif +#ifdef USE_OASIS + void initOasis(const bool writeOasisGrid); +#endif + private: TimePoint m_time; ConfigMap m_config; diff --git a/core/src/include/PrognosticData.hpp b/core/src/include/PrognosticData.hpp index 85426a0e8..2f1a5f06c 100644 --- a/core/src/include/PrognosticData.hpp +++ b/core/src/include/PrognosticData.hpp @@ -38,6 +38,8 @@ class PrognosticData : public ModelComponent, public Configured ModelState getState(const OutputLevel& lvl) const override { return getState(); } ModelState getStateRecursive(const OutputSpec& os) const override; + void setMetadata(const ModelMetadata& metadata); + static HelpMap& getHelpText(HelpMap& map, bool getAll); static HelpMap& getHelpRecursive(HelpMap& map, bool getAll); diff --git a/core/src/main.cpp b/core/src/main.cpp index bb94b4044..92f96c841 100644 --- a/core/src/main.cpp +++ b/core/src/main.cpp @@ -1,6 +1,6 @@ /*! * @file main.cpp - * @date 11 Aug 2021 + * @date 10 Sep 2024 * @author Tim Spain * @author Kacper Kornet */ @@ -9,6 +9,9 @@ #ifdef USE_MPI #include #endif +#ifdef USE_OASIS +#include +#endif #include "include/CommandLineParser.hpp" #include "include/ConfigurationHelpPrinter.hpp" @@ -21,6 +24,7 @@ int main(int argc, char* argv[]) { #ifdef USE_MPI MPI_Init(&argc, &argv); + MPI_Comm modelCommunicator = MPI_COMM_WORLD; #endif // USE_MPI // Pass the command line to Configurator to handle @@ -49,7 +53,15 @@ int main(int argc, char* argv[]) } else { // Construct the Model #ifdef USE_MPI - Nextsim::Model model(MPI_COMM_WORLD); +#ifdef USE_OASIS + /* We must call these oasis routines before any MPI communication takes place, to make sure + * we have the right communicator, i.e. modelCommunictor and not MPI_COMM_WORLD. */ + int compID; // Not actually used. Only useful for debugging + const std::string compName = "nextsim"; // Not useful for any setups we have in mind + OASIS_CHECK_ERR(oasis_c_init_comp(&compID, compName.c_str(), OASIS_COUPLED)); + OASIS_CHECK_ERR(oasis_c_get_localcomm(&modelCommunicator)); +#endif // USE_OASIS + Nextsim::Model model(modelCommunicator); #else Nextsim::Model model; #endif @@ -59,6 +71,9 @@ int main(int argc, char* argv[]) model.run(); } #ifdef USE_MPI +#ifdef USE_OASIS + OASIS_CHECK_ERR(oasis_c_terminate()); +#endif MPI_Finalize(); #endif diff --git a/core/test/ConfigOutput_test.cpp b/core/test/ConfigOutput_test.cpp index 36a65f820..e284e0c42 100644 --- a/core/test/ConfigOutput_test.cpp +++ b/core/test/ConfigOutput_test.cpp @@ -1,7 +1,7 @@ /*! * @file ConfigOutput_test.cpp * - * @date 11 May 2023 + * @date 10 Sep 2024 * @author Tim Spain */ @@ -17,7 +17,6 @@ #include "include/FileCallbackCloser.hpp" #include "include/IStructure.hpp" #include "include/ModelArray.hpp" -#include "include/ModelArrayRef.hpp" #include "include/ModelComponent.hpp" #include "include/ModelMetadata.hpp" #include "include/ModelState.hpp" @@ -30,8 +29,8 @@ #include #include -#include #include +#include const std::string test_files_dir = TEST_FILES_DIR; #ifdef USE_MPI @@ -54,10 +53,10 @@ TEST_CASE("Test periodic output") #ifdef USE_MPI if (test_rank == 0) { - ModelArray::setDimension(ModelArray::Dimension::X, nx, 1, 0); + ModelArray::setDimension(ModelArray::Dimension::X, nx, 1, 0); } if (test_rank == 1) { - ModelArray::setDimension(ModelArray::Dimension::X, nx, 1, 1); + ModelArray::setDimension(ModelArray::Dimension::X, nx, 1, 1); } ModelArray::setDimension(ModelArray::Dimension::Y, ny, ny, 0); ModelArray::setDimension(ModelArray::Dimension::Z, NZLevels::get(), NZLevels::get(), 0); diff --git a/physics/src/include/OASISCoupled.hpp b/physics/src/include/OASISCoupled.hpp new file mode 100644 index 000000000..a4cdf8191 --- /dev/null +++ b/physics/src/include/OASISCoupled.hpp @@ -0,0 +1,38 @@ +/*! + * @file OASISCoupled.hpp + * + * @date 09 Sep 2024 + * @author Einar Ólason + */ + +#ifndef OASISCOUPLED_HPP +#define OASISCOUPLED_HPP + +#ifdef USE_OASIS +#include +#endif + +namespace Nextsim { + +class OASISCoupled { +public: + virtual std::string getName() const { return "OASISCoupled"; } + +#ifdef USE_OASIS + int OASISTime; + + // Set the "OASIS time" (seconds since start) to zero + OASISCoupled() { OASISTime = 0; } + + // Increment the "OASIS" time by the number of seconds in the time step + // Could be any time unit + // Must be called at the end of the child class' update or updateAfter call. + void updateOASISTime(const TimestepTime& tst) { OASISTime += tst.step.seconds(); } +#else + const std::string OASISError = std::string(": OASIS support not compiled in.\n"); +#endif +}; + +} + +#endif // OASISCOUPLED_HPP diff --git a/physics/src/modules/OceanBoundaryModule/OASISCoupledOcean.cpp b/physics/src/modules/OceanBoundaryModule/OASISCoupledOcean.cpp new file mode 100644 index 000000000..8a580bab4 --- /dev/null +++ b/physics/src/modules/OceanBoundaryModule/OASISCoupledOcean.cpp @@ -0,0 +1,201 @@ +/*! + * @file OASISCoupledOcean.cpp + * + * @date 10 Sep 2024 + * @author Tim Spain + * @author Einar Ólason + */ + +#include "include/OASISCoupledOcean.hpp" +#include "include/IIceOceanHeatFlux.hpp" +#include "include/Module.hpp" +#include "include/constants.hpp" + +namespace Nextsim { + +void OASISCoupledOcean::setMetadata(const ModelMetadata& metadata) +{ +#ifdef USE_OASIS + // OASIS defining variable + + /* Populate the couplingId map with the id string and number pair. We need to do this seperately + * for the input (get) and output (put) variables. */ + // TODO: coplingID should be a map of where the pair is idNumber and + // pointer to the ModelArray + for (std::string idString : cplStringsIn) { + int idNumber; + OASIS_CHECK_ERR(oasis_c_def_var(&idNumber, idString.c_str(), metadata.OASISPartitionId, + bundleSize, OASIS_IN, OASIS_DOUBLE)); + couplingId[idString] = idNumber; + } + + for (std::string idString : cplStringsOut) { + int idNumber; + OASIS_CHECK_ERR(oasis_c_def_var(&idNumber, idString.c_str(), metadata.OASISPartitionId, + bundleSize, OASIS_OUT, OASIS_DOUBLE)); + couplingId[idString] = idNumber; + } +#else + throw std::runtime_error(std::string(__func__) + ": " + OASISError); +#endif +} + +void OASISCoupledOcean::updateBefore(const TimestepTime& tst) +{ + // Directly set the array values +#ifdef USE_OASIS + + int kinfo; + + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(SSTKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &sst[0], &kinfo)); + + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(SSSKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &sss[0], &kinfo)); + + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(UOceanKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &u[0], &kinfo)); + + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(VOceanKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &v[0], &kinfo)); + + // TODO: Implement ssh reading and passing to dynamics! + // OASIS_CHECK_ERR(oasis_c_get(couplingId.at(SSHKey), OASISTime, dimension0, dimension1, + // bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &ssh[0], &kinfo)); + + if (couplingId.find(MLDKey) != couplingId.end()) { + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(MLDKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &mld[0], &kinfo)); + } else { + mld = firstLayerDepth; + } + + cpml = Water::cp * Water::rho * mld; + + overElements( + std::bind(&OASISCoupledOcean::updateTf, this, std::placeholders::_1, std::placeholders::_2), + TimestepTime()); + + Module::getImplementation().update(tst); + +#else + throw std::runtime_error(std::string(__func__) + ": " + OASISError); +#endif +} + +void OASISCoupledOcean::updateAfter(const TimestepTime& tst) +{ +#ifdef USE_OASIS + int kinfo; + + // TODO We still need the the actual data + HField dummy; + dummy.resize(); + dummy.setData(184.); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(TauXKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(TauYKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(EMPKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(QSWKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(QNoSunKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(SFluxKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(TauModKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(CIceKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo)); + + // Increment the "OASIS" time by the number of seconds in the time step + updateOASISTime(tst); +#else + throw std::runtime_error(std::string(__func__) + ": " + OASISError); +#endif +} + +void OASISCoupledOcean::configure() +{ + SSTKey = Configured::getConfiguration(SSTConfigKey, SSTKeyDefault); + SSSKey = Configured::getConfiguration(SSSConfigKey, SSSKeyDefault); + UOceanKey = Configured::getConfiguration(UOceanConfigKey, UOceanKeyDefault); + VOceanKey = Configured::getConfiguration(VOceanConfigKey, VOceanKeyDefault); + SSHKey = Configured::getConfiguration(SSHConfigKey, SSHKeyDefault); + MLDKey = Configured::getConfiguration(MLDConfigKey, MLDKeyDefault); + TauXKey = Configured::getConfiguration(TauXConfigKey, TauXKeyDefault); + TauYKey = Configured::getConfiguration(TauYConfigKey, TauYKeyDefault); + TauModKey = Configured::getConfiguration(TauModConfigKey, TauModKeyDefault); + EMPKey = Configured::getConfiguration(EMPConfigKey, EMPKeyDefault); + QNoSunKey = Configured::getConfiguration(QNoSunConfigKey, QNoSunKeyDefault); + QSWKey = Configured::getConfiguration(QSWConfigKey, QSWKeyDefault); + SFluxKey = Configured::getConfiguration(SFluxConfigKey, SFluxKeyDefault); + CIceKey = Configured::getConfiguration(CIceConfigKey, CIceKeyDefault); + + firstLayerDepth = Configured::getConfiguration(layerDepthConfigKey, FIRST_LAYER_DEPTH); + + cplStringsIn + = { SSTKey, SSSKey, UOceanKey, VOceanKey, SSHKey }; + if (Configured::getConfiguration(exchangeFirstLayerConfigKey, EXCHANGE_FIRST_LAYER)) { + cplStringsIn.push_back(MLDKey); + } + cplStringsOut + = { TauXKey, TauYKey, TauModKey, EMPKey, QNoSunKey, QSWKey, SFluxKey, CIceKey }; +} + +OASISCoupledOcean::HelpMap& OASISCoupledOcean::getHelpText(HelpMap& map, bool getAll) +{ + map[moduleName] = { + { layerDepthConfigKey, ConfigType::NUMERIC, { "0", "∞" }, std::to_string(FIRST_LAYER_DEPTH), + "m", "Depth of the first ocean model layer (if this is fixed)." }, + { exchangeFirstLayerConfigKey, ConfigType::BOOLEAN, { "true", "false" }, + std::to_string(EXCHANGE_FIRST_LAYER), "", + "Use the thickness of the first ocean layer provided through the coupler" }, + { SSTConfigKey, ConfigType::STRING, {}, SSTKeyDefault, "", + "The field name for sea surface temperature used in namcouple" }, + { SSSConfigKey, ConfigType::STRING, {}, SSSKeyDefault, "", + "The field name for sea surface salinity used in namcouple" }, + { UOceanConfigKey, ConfigType::STRING, {}, UOceanKeyDefault, "", + "The field name for ocean u-velocity used in namcouple" }, + { VOceanConfigKey, ConfigType::STRING, {}, VOceanKeyDefault, "", + "The field name for ocean v-velocity used in namcouple" }, + { SSHConfigKey, ConfigType::STRING, {}, SSHKeyDefault, "", + "The field name for sea surface height used in namcouple" }, + { MLDConfigKey, ConfigType::STRING, {}, MLDKeyDefault, "", + "The field name for the thickness of the first ocean layer in namcouple (if that's " + "defined)." }, + { TauXConfigKey, ConfigType::STRING, {}, TauXKeyDefault, "", + "The field name for the x-component of ice-ocean stress in namcouple." }, + { TauYConfigKey, ConfigType::STRING, {}, TauYKeyDefault, "", + "The field name for the y-component of ice-ocean stress in namcouple." }, + { TauModConfigKey, ConfigType::STRING, {}, TauModKeyDefault, "", + "The field name for the modulus of ice-ocean stress in namcouple." }, + { EMPConfigKey, ConfigType::STRING, {}, EMPKeyDefault, "", + "The field name for freshwater flux used in namcouple" }, + { QNoSunConfigKey, ConfigType::STRING, {}, QNoSunKeyDefault, "", + "The field name for non-solar flux used in namcouple" }, + { QSWConfigKey, ConfigType::STRING, {}, QSWKeyDefault, "", + "The field name for showrt-wave flux used in namcouple" }, + { QSWConfigKey, ConfigType::STRING, {}, QSWKeyDefault, "", + "The field name for showrt-wave flux used in namcouple" }, + { SFluxConfigKey, ConfigType::STRING, {}, SFluxKeyDefault, "", + "The field name for salt flux used in namcouple" }, + { CIceConfigKey, ConfigType::STRING, {}, CIceKeyDefault, "", + "The field name for sea-ice concentration used in namcouple" }, + }; + return map; +} +OASISCoupledOcean::HelpMap& OASISCoupledOcean::getHelpRecursive(HelpMap& map, bool getAll) +{ + return getHelpText(map, getAll); +} + +} /* namespace Nextsim */ diff --git a/physics/src/modules/OceanBoundaryModule/include/OASISCoupledOcean.hpp b/physics/src/modules/OceanBoundaryModule/include/OASISCoupledOcean.hpp new file mode 100644 index 000000000..02cc37be6 --- /dev/null +++ b/physics/src/modules/OceanBoundaryModule/include/OASISCoupledOcean.hpp @@ -0,0 +1,106 @@ +/*! + * @file OASISCoupledOcean.hpp + * + * @date 09 Sep 2024 + * @author Tim Spain + * @author Einar Ólason + */ + +#ifndef OASISCOUPLEDOCEAN_HPP +#define OASISCOUPLEDOCEAN_HPP + +#include "include/IFreezingPoint.hpp" +#include "include/IOceanBoundary.hpp" +#include "include/Module.hpp" +#include "include/OASISCoupled.hpp" +#include "include/SlabOcean.hpp" + +namespace Nextsim { + +static const std::string moduleName = "OASISCoupledOcean"; + +static const std::string layerDepthConfigKey = moduleName + ".layer_depth"; +static const std::string exchangeFirstLayerConfigKey = moduleName + ".exchange_first_layer"; + +static const double FIRST_LAYER_DEPTH = 1.; // There really is no sensible default here(!) +static const bool EXCHANGE_FIRST_LAYER = false; + +static const std::string SSTKeyDefault = "I_SST"; +static const std::string SSSKeyDefault = "I_SSS"; +static const std::string UOceanKeyDefault = "I_Uocn"; +static const std::string VOceanKeyDefault = "I_Vocn"; +static const std::string SSHKeyDefault = "I_SSH"; +static const std::string MLDKeyDefault = "I_MLD"; // This one is optional +static const std::string TauXKeyDefault = "I_taux"; +static const std::string TauYKeyDefault = "I_tauy"; +static const std::string TauModKeyDefault = "I_taumod"; +static const std::string EMPKeyDefault = "I_fwflux"; +static const std::string QNoSunKeyDefault = "I_rsnos"; +static const std::string QSWKeyDefault = "I_rsso"; +static const std::string SFluxKeyDefault = "I_sfi"; +static const std::string CIceKeyDefault = "I_conc"; + +static const std::string SSTConfigKey = ".sea_surface_temperature"; +static const std::string SSSConfigKey = ".sea_surface_salinity"; +static const std::string UOceanConfigKey = ".ocean_u_velocity"; +static const std::string VOceanConfigKey = ".ocean_v_velocity"; +static const std::string SSHConfigKey = ".sea_surface_height"; +static const std::string MLDConfigKey = ".first_ocean_layer_depth"; // This one is optional +static const std::string TauXConfigKey = ".ice_ocean_stress_x"; +static const std::string TauYConfigKey = ".ice_ocean_stress_y"; +static const std::string TauModConfigKey = ".ice_ocean_stress_modulo"; +static const std::string EMPConfigKey = ".fresh_water_flux"; +static const std::string QNoSunConfigKey = ".non_solar_heatflux"; +static const std::string QSWConfigKey = ".short_wave_flux"; +static const std::string SFluxConfigKey = ".salt_flux"; +static const std::string CIceConfigKey = ".sea_ice_concentration"; + +//* Ocean boundary data values that are hardcoded. +class OASISCoupledOcean : public IOceanBoundary, + public OASISCoupled, + public Configured { +public: + OASISCoupledOcean() + : IOceanBoundary() + { + } + ~OASISCoupledOcean() { OASISCoupled::~OASISCoupled(); } + + std::string getName() const override { return moduleName; } + void updateBefore(const TimestepTime& tst) override; + void updateAfter(const TimestepTime& tst) override; + void setMetadata(const ModelMetadata& metadata) override; + + void configure() override; + + static HelpMap& getHelpText(HelpMap& map, bool getAll); + static HelpMap& getHelpRecursive(HelpMap& map, bool getAll); + +private: + int bundleSize = 1; // Always "unbundled", as per the OASIS manual + double firstLayerDepth = FIRST_LAYER_DEPTH; + + const int dimension0 = ModelArray::dimensions(ModelArray::Type::H)[0]; + const int dimension1 = ModelArray::dimensions(ModelArray::Type::H)[1]; + + SlabOcean slabOcean; + + void updateTf(size_t i, const TimestepTime& tst) + { + tf[i] = Module::getImplementation()(sss[i]); + } + + // A map to relate the strings in the namcouple file to the numbers def_var spits out + std::map couplingId; + std::string SSTKey, SSSKey, UOceanKey, VOceanKey, SSHKey, MLDKey, TauXKey, TauYKey, TauModKey, + EMPKey, QNoSunKey, QSWKey, SFluxKey, CIceKey; + + std::vector cplStringsIn + = { SSTKey, SSSKey, UOceanKey, VOceanKey, SSHKey }; // MLDKey can be added to this one + std::vector cplStringsOut + = { TauXKey, TauYKey, TauModKey, EMPKey, QNoSunKey, QSWKey, SFluxKey, CIceKey }; +}; + +} /* namespace Nextsim */ + +#endif /* OASISCOUPLEDOCEAN_HPP */ diff --git a/physics/src/modules/OceanBoundaryModule/module.cfg b/physics/src/modules/OceanBoundaryModule/module.cfg index ddadaf128..312dc572a 100644 --- a/physics/src/modules/OceanBoundaryModule/module.cfg +++ b/physics/src/modules/OceanBoundaryModule/module.cfg @@ -31,3 +31,8 @@ has_help = true file_prefix = BenchmarkOcean description = Calculated ocean for the DG dynamics benchmark. has_help = false + +[Nextsim::OASISCoupledOcean] +file_prefix = OASISCoupledOcean +description = Receive and send ocean fields via the OASIS coupler +has_help = false diff --git a/physics/src/modules/include/IAtmosphereBoundary.hpp b/physics/src/modules/include/IAtmosphereBoundary.hpp index 5b64a734a..e62a4a80c 100644 --- a/physics/src/modules/include/IAtmosphereBoundary.hpp +++ b/physics/src/modules/include/IAtmosphereBoundary.hpp @@ -7,6 +7,7 @@ #include "include/ModelArrayRef.hpp" #include "include/ModelComponent.hpp" +#include "include/ModelMetadata.hpp" #include "include/Time.hpp" #ifndef IATMOSPHEREBOUNDARY_HPP @@ -15,12 +16,12 @@ namespace Nextsim { namespace CouplingFields { -constexpr TextTag SUBL = "SUBL"; // sublimation mass flux kg s⁻¹ m⁻² -constexpr TextTag SNOW = "SNOW"; // snowfall mass flux kg s⁻¹ m⁻² -constexpr TextTag RAIN = "RAIN"; // rainfall mass flux kg s⁻¹ m⁻² -constexpr TextTag EVAP = "EVAP"; // evaporation mass flux kg s⁻¹ m⁻² -constexpr TextTag WIND_U = "WIND_U"; // x-aligned wind component m s⁻¹ -constexpr TextTag WIND_V = "WIND_V"; // y-aligned wind component m s⁻¹ + constexpr TextTag SUBL = "SUBL"; // sublimation mass flux kg s⁻¹ m⁻² + constexpr TextTag SNOW = "SNOW"; // snowfall mass flux kg s⁻¹ m⁻² + constexpr TextTag RAIN = "RAIN"; // rainfall mass flux kg s⁻¹ m⁻² + constexpr TextTag EVAP = "EVAP"; // evaporation mass flux kg s⁻¹ m⁻² + constexpr TextTag WIND_U = "WIND_U"; // x-aligned wind component m s⁻¹ + constexpr TextTag WIND_V = "WIND_V"; // y-aligned wind component m s⁻¹ } //! An interface class for the atmospheric inputs into the ice physics. @@ -29,15 +30,15 @@ class IAtmosphereBoundary : public ModelComponent { IAtmosphereBoundary() : qia(ModelArray::Type::H) , dqia_dt(ModelArray::Type::H) - , qow(ModelArray::Type::H) - , subl(ModelArray::Type::H) - , snow(ModelArray::Type::H) - , rain(ModelArray::Type::H) - , evap(ModelArray::Type::H) - , emp(ModelArray::Type::H) - , uwind(ModelArray::Type::U) - , vwind(ModelArray::Type::V) - , penSW(ModelArray::Type::H) + , qow(ModelArray::Type::H) + , subl(ModelArray::Type::H) + , snow(ModelArray::Type::H) + , rain(ModelArray::Type::H) + , evap(ModelArray::Type::H) + , emp(ModelArray::Type::H) + , uwind(ModelArray::Type::U) + , vwind(ModelArray::Type::V) + , penSW(ModelArray::Type::H) { m_couplingArrays.registerArray(CouplingFields::SUBL, &subl, RW); m_couplingArrays.registerArray(CouplingFields::SNOW, &snow, RW); @@ -62,6 +63,7 @@ class IAtmosphereBoundary : public ModelComponent { ModelState getState(const OutputLevel&) const override { return getState(); } std::string getName() const override { return "IAtmosphereBoundary"; } + virtual void setMetadata(const ModelMetadata& metadata) { } void setData(const ModelState::DataMap& ms) override { qia.resize(); @@ -79,7 +81,6 @@ class IAtmosphereBoundary : public ModelComponent { virtual void update(const TimestepTime& tst) { } protected: - ModelArrayReferenceStore& couplingArrays() { return m_couplingArrays; } HField qia; diff --git a/physics/src/modules/include/IOceanBoundary.hpp b/physics/src/modules/include/IOceanBoundary.hpp index b85c4cb06..fa1a942dc 100644 --- a/physics/src/modules/include/IOceanBoundary.hpp +++ b/physics/src/modules/include/IOceanBoundary.hpp @@ -9,15 +9,16 @@ #define IOCEANBOUNDARY_HPP #include "include/ModelComponent.hpp" +#include "include/ModelMetadata.hpp" namespace Nextsim { namespace CouplingFields { -constexpr TextTag SST = "SST"; // sea surface temperature ˚C -constexpr TextTag SSS = "SSS"; // sea surface salinity PSU -constexpr TextTag MLD = "MLD"; // Mixed layer or slab ocean depth m -constexpr TextTag OCEAN_U = "U"; // x(east)-ward ocean current m s⁻¹ -constexpr TextTag OCEAN_V = "V"; // y(north)-ward ocean current m s⁻¹ + constexpr TextTag SST = "SST"; // sea surface temperature ˚C + constexpr TextTag SSS = "SSS"; // sea surface salinity PSU + constexpr TextTag MLD = "MLD"; // Mixed layer or slab ocean depth m + constexpr TextTag OCEAN_U = "U"; // x(east)-ward ocean current m s⁻¹ + constexpr TextTag OCEAN_V = "V"; // y(north)-ward ocean current m s⁻¹ } //! An interface class for the oceanic inputs into the ice physics. class IOceanBoundary : public ModelComponent { @@ -44,6 +45,7 @@ class IOceanBoundary : public ModelComponent { ModelState getState(const OutputLevel&) const override { return getState(); } std::string getName() const override { return "IOceanBoundary"; } + virtual void setMetadata(const ModelMetadata& metadata) { } void setData(const ModelState::DataMap& ms) override { qio.resize(); diff --git a/physics/test/CMakeLists.txt b/physics/test/CMakeLists.txt index 2b01dcc4d..0a204ad34 100644 --- a/physics/test/CMakeLists.txt +++ b/physics/test/CMakeLists.txt @@ -37,6 +37,12 @@ if(ENABLE_MPI) PRIVATE USE_MPI TEST_FILES_DIR=\"${CMAKE_CURRENT_SOURCE_DIR}\" ) target_link_libraries(testTOPAZOcn_MPI1 PRIVATE nextsimlib doctest::doctest) + if (ENABLE_OASIS) + add_executable(testOASISCoupledOcean_MPI2 "OASISCoupledOcean_test.cpp" "MainMPI.cpp") + target_include_directories(testOASISCoupledOcean_MPI2 PRIVATE "${ModulesRoot}/OceanBoundaryModule") + target_compile_definitions(testOASISCoupledOcean_MPI2 PRIVATE TEST_FILES_DIR=\"${CMAKE_CURRENT_SOURCE_DIR}\") + target_link_libraries(testOASISCoupledOcean_MPI2 PRIVATE nextsimlib doctest::doctest) + endif () else() add_executable(testERA5Atm "ERA5Atm_test.cpp") target_compile_definitions(testERA5Atm PRIVATE TEST_FILES_DIR=\"${CMAKE_CURRENT_SOURCE_DIR}\") diff --git a/physics/test/MainMPI.cpp b/physics/test/MainMPI.cpp new file mode 100644 index 000000000..5fde9684e --- /dev/null +++ b/physics/test/MainMPI.cpp @@ -0,0 +1,30 @@ +/*! + * @file MainMPI.cpp + * + * @date 13 Sep 2024 + * @author Einar Ólason + * + * This file is required so that the MPI enabled doc-test OASISCoupledOcean_test.cpp (and any + * others) can run. + */ + +#define DOCTEST_CONFIG_IMPLEMENT + +#include + +int main(int argc, char** argv) +{ + doctest::mpi_init_thread(argc, argv, MPI_THREAD_MULTIPLE); + + doctest::Context ctx; + ctx.setOption("reporters", "MpiConsoleReporter"); + ctx.setOption("reporters", "MpiFileReporter"); + ctx.setOption("force-colors", true); + ctx.applyCommandLine(argc, argv); + + int test_result = ctx.run(); + + doctest::mpi_finalize(); + + return test_result; +} diff --git a/physics/test/OASISCoupledOcean_test.cpp b/physics/test/OASISCoupledOcean_test.cpp new file mode 100644 index 000000000..7a8f32d30 --- /dev/null +++ b/physics/test/OASISCoupledOcean_test.cpp @@ -0,0 +1,66 @@ +/*! + * @file OASISCoupledOcean_test.cpp + * + * @date 13 Sep 2024 + * @author Einar Ólason + */ + +#include + +#include "include/OASISCoupledOcean.hpp" + +namespace Nextsim { + +TEST_SUITE_BEGIN("OASISCoupledOcean"); +MPI_TEST_CASE("OASIS init put and get", 1) +{ + MPI_Comm modelCommunicator; + int compID; // Not actually used. Only useful for debugging + const std::string compName = "nextsim"; // Not useful for any setups we have in mind + OASIS_CHECK_ERR(oasis_c_init_comp(&compID, compName.c_str(), OASIS_COUPLED)); + OASIS_CHECK_ERR(oasis_c_get_localcomm(&modelCommunicator)); + + ModelArray::setDimensions(ModelArray::Type::H, { 1, 1 }); + ModelArray::setDimensions(ModelArray::Type::Z, { 1, 1, 1 }); + + double sstIn = -1.84; + double sssIn = 28.0; + //double mldIn = 14.8; + double uIn = -0.14; + double vIn = 0.71; + + HField cice(ModelArray::Type::H); + cice = 1.0; + ModelComponent::getStore().registerArray(Protected::C_ICE, &cice, RO); + OASISCoupledOcean ocpl; + ModelMetadata metadata; + const std::vector partInfo = { OASIS_Serial, 1, 1 }; + OASIS_CHECK_ERR(oasis_c_def_partition( + &metadata.OASISPartitionId, OASIS_Serial_Params, &partInfo[0], OASIS_No_Gsize, compName.c_str())); + + ocpl.setData(ModelState::DataMap()); + ocpl.configure(); + ocpl.setMetadata(metadata); + OASIS_CHECK_ERR(oasis_c_enddef()); + + ocpl.updateBefore(TimestepTime()); + ModelArrayRef sst(ModelComponent::getStore()); + ModelArrayRef sss(ModelComponent::getStore()); + ModelArrayRef u(ModelComponent::getStore()); + ModelArrayRef v(ModelComponent::getStore()); + std::cout << "Received SST at time " << ocpl.OASISTime << ": " << sst[0] << std::endl ; + std::cout << "Received SSS at time " << ocpl.OASISTime << ": " << sss[0] << std::endl ; + std::cout << "Received OCEAN_U at time " << ocpl.OASISTime << ": " << u[0] << std::endl ; + std::cout << "Received OCEAN_V at time " << ocpl.OASISTime << ": " << v[0] << std::endl ; + REQUIRE(sst[0] == sstIn); + REQUIRE(sss[0] == sssIn); + REQUIRE(u[0] == uIn); + REQUIRE(v[0] == vIn); + //REQUIRE(mld[0] == mldIn); + + ocpl.updateAfter(TimestepTime()); + + OASIS_CHECK_ERR(oasis_c_terminate()); +} +TEST_SUITE_END(); +} diff --git a/physics/test/generate_OASIS_NetCDF.sh b/physics/test/generate_OASIS_NetCDF.sh new file mode 100755 index 000000000..a7b265826 --- /dev/null +++ b/physics/test/generate_OASIS_NetCDF.sh @@ -0,0 +1,75 @@ +#!/usr/bin/bash + +InVars=(I_SST:-1.84 I_SSS:28 I_Uocn:-0.14 I_Vocn:0.71) +# Optionally add +# InVars+=(I_SSH:14.8 I_MLD:14.8) + +OutVars=(I_taux I_tauy I_taumod I_fwflux I_rsnos I_rsso I_sfi I_conc) + +cat > namcouple < tmp.cdl <> namcouple + +done + +cat >> namcouple <> namcouple <> namcouple