diff --git a/CMakeLists.txt b/CMakeLists.txt index 1f3a295f9..8474e92f1 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" @@ -155,6 +156,20 @@ if(ENABLE_XIOS) target_link_libraries(nextsimlib PUBLIC "${FORTRAN_RUNTIME_LIB}") 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) + target_compile_definitions(nextsimlib PUBLIC USE_OASIS) + + target_link_libraries(nextsimlib PUBLIC ${OASIS_LIBRARIES}) + 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/Dockerfiles/Dockerfile.devenv b/Dockerfiles/Dockerfile.devenv index b11682df4..cc9c7d721 100644 --- a/Dockerfiles/Dockerfile.devenv +++ b/Dockerfiles/Dockerfile.devenv @@ -50,6 +50,8 @@ COPY --from=builder /opt/views /opt/views # Create entrypoint script RUN { \ echo '#!/bin/sh'; \ + echo 'export OMPI_ALLOW_RUN_AS_ROOT=1'; \ + echo 'export OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1'; \ echo '. /opt/spack-environment/activate.sh'; \ echo 'exec "$@"'; \ } > /entrypoint.sh \ diff --git a/Dockerfiles/Dockerfile.oasis b/Dockerfiles/Dockerfile.oasis new file mode 100644 index 000000000..74d57a403 --- /dev/null +++ b/Dockerfiles/Dockerfile.oasis @@ -0,0 +1,30 @@ +FROM nextsimdg-dev-env:latest + +# A temproary Dockerfile with commands that should find their way into Dockerfile.devenv in the end + +# Clone the repository +RUN git clone https://gitlab.com/cerfacs/oasis3-mct.git +WORKDIR /oasis3-mct/util/make_dir + +# Set some ENV's and copy the make file include +ENV OASIS_COUPLE=/oasis3-mct +ENV OASIS_ENV=oasis_docker_spack_gcc +ENV OASIS_DIR=$OASIS_COUPLE/build.$OASIS_ENV +COPY make.oasis_docker make.$OASIS_ENV + +# make! +RUN make -f TopMakefileOasis3 static-libs shared-libs pyoasis + +# Install python packages that pyoasis may need +RUN apt-get install -y python3-mpi4py python3-isodate + +# Append the OASIS library directory to LD_LIBRARY_PATH and append it to PYTHONPATH +# We assume the last line of the file is 'exec "$@"' and try to retain that +RUN sed -i~ '$d' /entrypoint.sh +RUN { \ + echo 'export LD_LIBRARY_PATH=$LD_LIBRARY_PATH':$OASIS_DIR/lib; \ + echo 'export PYTHONPATH=$PYTHONPATH':$OASIS_DIR/python; \ + echo 'exec "$@"'; \ + } >> /entrypoint.sh + +WORKDIR / \ No newline at end of file diff --git a/Dockerfiles/make.oasis_docker b/Dockerfiles/make.oasis_docker new file mode 100644 index 000000000..c611b07cf --- /dev/null +++ b/Dockerfiles/make.oasis_docker @@ -0,0 +1,82 @@ +# +# Include file for OASIS3 Makefile for a Linux docker image using spack and gcc +# +############################################################################### +# +# CHAN : communication technique used in OASIS3 (MPI1/MPI2) +CHAN = MPI1 +# +# Paths for libraries, object files and binaries +# +# COUPLE : path for oasis3-mct main directory +COUPLE := $(OASIS_COUPLE) +# +# ARCHDIR : directory created when compiling +ARCHDIR := $(COUPLE)/build.$(OASIS_ENV) +# +# MPI library +MPIDIR = /opt/view/ +MPIBIN = $(MPIDIR)/bin +MPI_INCLUDE = $(MPIDIR)/include +MPILIB = -L$(MPIDIR)/lib -lmpi +MPIRUN = $(MPIBIN)/mpirun +# +# NETCDF library of the system +NETCDF_INCLUDE = /opt/view/include +NETCDF_LIBRARY = -L/opt/view/lib -Wl,-rpath,/opt/view/lib -lnetcdff -lnetcdf -lm +# +# Compiling and other commands +MAKE = gmake +F90 = $(MPIBIN)/mpif90 -I$(MPI_INCLUDE) +F = $(F90) +f90 = $(F90) +f = $(F90) +CC = $(MPIBIN)/mpicc -I$(MPI_INCLUDE) +LD = $(MPIBIN)/mpif90 $(MPILIB) +DYNOPT = -fPIC +LDDYNOPT = -shared +AR = ar +ARFLAGS = -ruv +# Fortran libraries for C linking +F2C_LIBS = -lmpi_mpifh -lgfortran -lm +# +# CPP keys and compiler options +# +CPPDEF = -Duse_comm_$(CHAN) -D__VERBOSE -DTREAT_OVERLAY +# +FCBASEFLAGS := -fallow-argument-mismatch -ffree-line-length-0 +ifeq ($(OASIS_DEBUG), ) + ifeq ($(OASIS_OPENMP), ) + FCBASEFLAGS += -O2 + else + FCBASEFLAGS += -O2 -fopenmp + endif +else + ifeq ($(OASIS_OPENMP), ) + FCBASEFLAGS += -g -fbounds-check + else + FCBASEFLAGS += -g -fbounds-check -fopenmp + endif +endif +CCBASEFLAGS := +ifeq ($(OASIS_OPENMP), ) + CCBASEFLAGS += +else + CCBASEFLAGS += -fopenmp +endif +# +# INC_DIR : includes all *mod for each library + INC_DIR = -I$(ARCHDIR)/include +# FLIBS : for toys when linking in local Makefile + FLIBS=${NETCDF_LIBRARY} +################### +# +F90FLAGS = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE) +f90FLAGS = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE) +FFLAGS = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE) +fFLAGS = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE) +CCFLAGS = $(CCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE) +LDFLAGS = $(FCBASEFLAGS) +F2C_LDFLAGS = $(F2C_LIBS) +# +############################################################################# diff --git a/cmake/FindOASIS.cmake b/cmake/FindOASIS.cmake new file mode 100644 index 000000000..19be22304 --- /dev/null +++ b/cmake/FindOASIS.cmake @@ -0,0 +1,31 @@ +# Find oasis +# +# Please pass the -DOASIS_DIR variable to cmake (location of the OASIS libraries). + +# Search for all library files individually and merge to a single OASIS_LIBRARIES variable for the parent CMakeLists.txt +find_library(LIB_OASIS NAMES oasis.cbind HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH) +find_library(LIB_MCT NAMES mct HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH) +find_library(LIB_MPEU NAMES mpeu HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH) +find_library(LIB_PSMILE NAMES psmile.MPI1 psmile.MPI2 HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH) +find_library(LIB_SCRIP NAMES scrip HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH) + +# Only continue if all libraries are found +if (LIB_OASIS AND LIB_MCT AND LIB_MPEU AND LIB_PSMILE AND LIB_SCRIP) + set(OASIS_LIBRARIES ${LIB_OASIS} ${LIB_MCT} ${LIB_MPEU} ${LIB_PSMILE} ${LIB_SCRIP}) +endif () + +# Set OASIS_DIR (not strictly necessary, but nice to have) +get_filename_component(OASIS_LIB_DIR "${LIB_OASIS}" PATH) +set(OASIS_DIR "${OASIS_LIB_DIR}/../") +cmake_path(NORMAL_PATH OASIS_DIR) + +# The oasis include file should be where oasis.cbind lives +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_DIR 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/DevStep.cpp b/core/src/DevStep.cpp index 04f692fb7..f7dc28d6e 100644 --- a/core/src/DevStep.cpp +++ b/core/src/DevStep.cpp @@ -41,6 +41,8 @@ void DevStep::iterate(const TimestepTime& tst) mData.incrementTime(tst.step); if ((m_restartPeriod.seconds() > 0) && (mData.time() >= lastOutput + m_restartPeriod)) { std::string currentFileName = mData.time().format(m_restartFileName); + // Some MPI-IO implementations does not like colon in file names + std::replace(currentFileName.begin(), currentFileName.end(), ':', '_'); pData->writeRestartFile(currentFileName); lastOutput = mData.time(); } diff --git a/core/src/Model.cpp b/core/src/Model.cpp index 62012a242..74b35bc7c 100644 --- a/core/src/Model.cpp +++ b/core/src/Model.cpp @@ -40,6 +40,9 @@ static const std::map keyMap = { { 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 }; Model::Model() @@ -122,9 +125,15 @@ void Model::configure() // Get the coordinates from the ModelState for persistence metadata.extractCoordinates(initialState); +#ifdef USE_OASIS + const bool writeOasisGrid = Configured::getConfiguration(keyMap.at(WRITEOASISGRID_KEY), false); + metadata.initOasis(writeOasisGrid); +#endif + modelStep.setData(pData); modelStep.setRestartDetails(metadata.restartPeriod, metadata.finalFileName); pData.setData(initialState.data); + pData.setMetadata(metadata); } Model::HelpMap& Model::getHelpText(HelpMap& map, bool getAll) @@ -195,6 +204,8 @@ void Model::writeRestartFile() { auto& metadata = ModelMetadata::getInstance(); std::string formattedFileName = metadata.time().format(metadata.finalFileName); + // Some MPI-IO implementations does not like colon in file names + std::replace(formattedFileName.begin(), formattedFileName.end(), ':', '_'); pData.writeRestartFile(formattedFileName); } diff --git a/core/src/ModelArray.cpp b/core/src/ModelArray.cpp index ac2b8eb47..23ff6f3e9 100644 --- a/core/src/ModelArray.cpp +++ b/core/src/ModelArray.cpp @@ -175,6 +175,27 @@ ModelArray& ModelArray::clampBelow(const ModelArray& minArr) return *this; } +ModelArray ModelArray::sqrt() +{ + ModelArray sqrted = ModelArray(type); + sqrted.m_data.array() = m_data.array().sqrt(); + return sqrted; +} + +ModelArray ModelArray::sin() +{ + auto sined = ModelArray(type); + sined.m_data.array() = m_data.array().sin(); + return sined; +} + +ModelArray ModelArray::cos() +{ + auto cosed = ModelArray(type); + cosed.m_data.array() = m_data.array().cos(); + return cosed; +} + void ModelArray::setData(double value) { resize(); diff --git a/core/src/ModelMetadata.cpp b/core/src/ModelMetadata.cpp index 79c981c0b..2e3d7e4dc 100644 --- a/core/src/ModelMetadata.cpp +++ b/core/src/ModelMetadata.cpp @@ -25,6 +25,9 @@ #include #include #endif +#ifdef USE_OASIS +#include +#endif namespace Nextsim { @@ -191,6 +194,71 @@ 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 + void ModelMetadata::setTimes(const TimePoint& start, const TimePoint& stop, const Duration& step) { this->start = start; diff --git a/core/src/PrognosticData.cpp b/core/src/PrognosticData.cpp index 140d6d7cc..e0a25d57f 100644 --- a/core/src/PrognosticData.cpp +++ b/core/src/PrognosticData.cpp @@ -11,6 +11,10 @@ #include "include/NextsimModule.hpp" #include "include/gridNames.hpp" +#ifdef USE_OASIS +#include +#endif + namespace Nextsim { static const std::string checkFieldsKey = "debug.check_fields"; @@ -111,6 +115,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) { // Prepare everything diff --git a/core/src/include/Model.hpp b/core/src/include/Model.hpp index 79796732d..3b1ab44ea 100644 --- a/core/src/include/Model.hpp +++ b/core/src/include/Model.hpp @@ -43,6 +43,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 { return ConfigMap(); }; diff --git a/core/src/include/ModelArray.hpp b/core/src/include/ModelArray.hpp index 5f279a453..616661571 100644 --- a/core/src/include/ModelArray.hpp +++ b/core/src/include/ModelArray.hpp @@ -297,6 +297,18 @@ class ModelArray { * @param minArr the array of clamp minimum target values. */ ModelArray& clampBelow(const ModelArray& minArr); + /*! + * @brief Returns the per-element square root of the array + */ + ModelArray sqrt(); + /*! + * @brief Returns the per-element sin of the array + */ + ModelArray sin(); + /*! + * @brief Returns the per-element cos of the array + */ + ModelArray cos(); using MultiDim = std::vector; diff --git a/core/src/include/ModelMetadata.hpp b/core/src/include/ModelMetadata.hpp index 6952cd565..da190c26d 100644 --- a/core/src/include/ModelMetadata.hpp +++ b/core/src/include/ModelMetadata.hpp @@ -71,6 +71,10 @@ class ModelMetadata { } #endif +#ifdef USE_OASIS + int OASISPartitionId; +#endif + // finalize ModelMetadata static void finalize(); @@ -198,6 +202,10 @@ class ModelMetadata { neighbourArray neighbourHaloRecvPeriodic; #endif +#ifdef USE_OASIS + void initOasis(const bool writeOasisGrid); +#endif + std::string initialFileName; std::string finalFileName; // Period between restart file outputs diff --git a/core/src/include/PrognosticData.hpp b/core/src/include/PrognosticData.hpp index f4a2651e8..65d907baf 100644 --- a/core/src/include/PrognosticData.hpp +++ b/core/src/include/PrognosticData.hpp @@ -34,6 +34,8 @@ class PrognosticData : public CheckingModelComponent, public Configured #endif +#ifdef USE_OASIS +#include +#endif #include "include/CommandLineParser.hpp" #include "include/ConfigurationHelpPrinter.hpp" @@ -48,11 +51,20 @@ int main(int argc, char* argv[]) } else { // Construct the Model #ifdef USE_MPI - Nextsim::ModelMPI& modelMPI = Nextsim::ModelMPI::getInstance(MPI_COMM_WORLD); - Nextsim::Model model; +#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)); + MPI_Comm modelCommunicator; + OASIS_CHECK_ERR(oasis_c_get_localcomm(&modelCommunicator)); + Nextsim::ModelMPI& modelMPI = Nextsim::ModelMPI::getInstance(modelCommunicator); #else - Nextsim::Model model; + Nextsim::ModelMPI& modelMPI = Nextsim::ModelMPI::getInstance(MPI_COMM_WORLD); +#endif // USE_OASIS #endif + Nextsim::Model model; // Apply the model configuration model.configure(); // Run the Model @@ -62,6 +74,11 @@ int main(int argc, char* argv[]) return_code = -1; Nextsim::Logged::error(e.what()); } +#ifdef USE_MPI +#ifdef USE_OASIS + OASIS_CHECK_ERR(oasis_c_terminate()); +#endif +#endif } #ifdef USE_MPI MPI_Finalize(); diff --git a/core/src/modules/include/IDiagnosticOutput.hpp b/core/src/modules/include/IDiagnosticOutput.hpp index cdfb68d5c..5818dbbdd 100644 --- a/core/src/modules/include/IDiagnosticOutput.hpp +++ b/core/src/modules/include/IDiagnosticOutput.hpp @@ -23,7 +23,7 @@ class IDiagnosticOutput : public ModelComponent { */ #include "include/ProtectedArrayNames.ipp" #include "include/SharedArrayNames.ipp" - }) + }) { } virtual ~IDiagnosticOutput() = default; diff --git a/core/src/modules/include/IDynamics.hpp b/core/src/modules/include/IDynamics.hpp index 2be8c7ffd..86d76d4e7 100644 --- a/core/src/modules/include/IDynamics.hpp +++ b/core/src/modules/include/IDynamics.hpp @@ -86,6 +86,8 @@ class IDynamics : public ModelComponent { { uice.resize(); vice.resize(); + taux.resize(); + tauy.resize(); shear.resize(); divergence.resize(); sigmaI.resize(); diff --git a/core/test/ConfigOutput_test.cpp b/core/test/ConfigOutput_test.cpp index dbb2bb6e3..0130bc1bf 100644 --- a/core/test/ConfigOutput_test.cpp +++ b/core/test/ConfigOutput_test.cpp @@ -16,7 +16,6 @@ #include "include/Finalizer.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" diff --git a/core/test/ModelArray_test.cpp b/core/test/ModelArray_test.cpp index 52930d6e8..f01657f35 100644 --- a/core/test/ModelArray_test.cpp +++ b/core/test/ModelArray_test.cpp @@ -179,6 +179,9 @@ TEST_CASE("Arithmetic tests") OneDField negative = -rhs; REQUIRE(negative[0] == -3); REQUIRE(negative[1] == 5); + OneDField squareRoot = lhs.sqrt(); + REQUIRE(squareRoot[0] == 3.); + REQUIRE(squareRoot[1] == std::sqrt(lhs[1])); double three = 3; double four = 4; diff --git a/physics/src/include/OASISCoupled.hpp b/physics/src/include/OASISCoupled.hpp new file mode 100644 index 000000000..96876500e --- /dev/null +++ b/physics/src/include/OASISCoupled.hpp @@ -0,0 +1,66 @@ +/*! + * @file OASISCoupled.hpp + * + * @date 09 Sep 2024 + * @author Einar Ólason + */ + +#ifndef OASISCOUPLED_HPP +#define OASISCOUPLED_HPP + +#ifdef USE_OASIS +#include +#endif + +#include "include/ModelArray.hpp" +#include "include/Time.hpp" + +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 NoOASISError = std::string(": OASIS support not compiled in.\n"); +#endif + +protected: + void rotateVectorFromGreenland( + const HField& uIn, const HField& vIn, HField& uOut, HField& vOut) const + { + // Make sure the outputs have the right size + uOut.resize(); + vOut.resize(); + + uOut = uIn * cosAngle - vIn * sinAngle; + vOut = uIn * sinAngle + vIn * cosAngle; + }; + + void rotateVectorToGreenland( + const HField& uIn, const HField& vIn, HField& uOut, HField& vOut) const + { + // Make sure the outputs have the right size + uOut.resize(); + vOut.resize(); + + uOut = uIn * cosAngle + vIn * sinAngle; + vOut = -uIn * sinAngle + vIn * cosAngle; + }; + + HField cosAngle, sinAngle; +}; + +} + +#endif // OASISCOUPLED_HPP diff --git a/physics/src/include/SlabOcean.hpp b/physics/src/include/SlabOcean.hpp index 10d2930a2..0769ecaad 100644 --- a/physics/src/include/SlabOcean.hpp +++ b/physics/src/include/SlabOcean.hpp @@ -20,6 +20,8 @@ namespace Nextsim { */ class SlabOcean : public ModelComponent, public Configured { public: + SlabOcean() = delete; + SlabOcean(ModelArrayReferenceStore& coupingArrays) : qdw(ModelArray::Type::H) , fdw(ModelArray::Type::H) diff --git a/physics/src/modules/OceanBoundaryModule/OASISCoupledOcean.cpp b/physics/src/modules/OceanBoundaryModule/OASISCoupledOcean.cpp new file mode 100644 index 000000000..ba51cd329 --- /dev/null +++ b/physics/src/modules/OceanBoundaryModule/OASISCoupledOcean.cpp @@ -0,0 +1,269 @@ +/*! + * @file OASISCoupledOcean.cpp + * + * @author Tim Spain + * @author Einar Ólason + */ + +#include "include/OASISCoupledOcean.hpp" + +#include "include/Finalizer.hpp" +#include "include/IFreezingPoint.hpp" +#include "include/IIceOceanHeatFlux.hpp" +#include "include/NextsimModule.hpp" + +namespace Nextsim { + +void OASISCoupledOcean::setData(const ModelState::DataMap& ms) +{ + IOceanBoundary::setData(ms); + + // Precalculate cos and sin of the azimut angle to use in the rotateVector... functions + HField gridAzimuthAngle = ms.at(gridAzimuthName); + cosAngle = gridAzimuthAngle.cos(); + sinAngle = gridAzimuthAngle.sin(); +} + +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 separately + * 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; + + Logged::debug("OASISCoupledOcean: " + idString + " has id " + std::to_string(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; + + Logged::debug("OASISCoupledOcean: " + idString + " has id " + std::to_string(idNumber)); + } +#else + throw std::runtime_error(std::string(__func__) + ": " + NoOASISError); +#endif +} + +void OASISCoupledOcean::updateBefore(const TimestepTime& tst) +{ + // Directly set the array values +#ifdef USE_OASIS + + int kinfo; + const int dimension0 = ModelArray::dimensions(ModelArray::Type::H)[0]; + const int dimension1 = ModelArray::dimensions(ModelArray::Type::H)[1]; + + Logged::debug("OASISCoupledOcean: Receiving " + SSTKey + " with id " + + std::to_string(couplingId.at(SSTKey))); + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(SSTKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &sst[0], &kinfo)); + + Logged::debug("OASISCoupledOcean: Receiving " + SSSKey + " with id " + + std::to_string(couplingId.at(SSSKey))); + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(SSSKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &sss[0], &kinfo)); + + HField uOnGrid; + uOnGrid.resize(); + Logged::debug("OASISCoupledOcean: Receiving " + UOceanKey + " with id " + + std::to_string(couplingId.at(UOceanKey))); + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(UOceanKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &uOnGrid[0], &kinfo)); + + HField vOnGrid; + vOnGrid.resize(); + Logged::debug("OASISCoupledOcean: Receiving " + VOceanKey + " with id " + + std::to_string(couplingId.at(VOceanKey))); + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(VOceanKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &vOnGrid[0], &kinfo)); + + Logged::debug("OASISCoupledOcean: Receiving " + SSHKey + " with id " + + std::to_string(couplingId.at(SSHKey))); + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(SSHKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &ssh[0], &kinfo)); + + Logged::debug("OASISCoupledOcean: Receiving " + FracQSWKey + " with id " + + std::to_string(couplingId.at(FracQSWKey))); + OASIS_CHECK_ERR(oasis_c_get(couplingId.at(FracQSWKey), OASISTime, dimension0, dimension1, + bundleSize, OASIS_DOUBLE, OASIS_COL_MAJOR, &fracQSWAbs[0], &kinfo)); + + if (couplingId.find(MLDKey) != couplingId.end()) { + Logged::debug("OASISCoupledOcean: Receiving " + MLDKey + " with id " + + std::to_string(couplingId.at(MLDKey))); + 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::rhoOcean * Water::cp * mld; + rotateVectorToGreenland(uOnGrid, vOnGrid, u, v); + + overElements([this](size_t i, const TimestepTime& tsTime) { this->updateTf(i, tsTime); }, tst); + + Module::getImplementation().update(tst); + +#else + throw std::runtime_error(std::string(__func__) + ": " + NoOASISError); +#endif +} + +void OASISCoupledOcean::updateAfter(const TimestepTime& tst) +{ +#ifdef USE_OASIS + mergeFluxes(tst); + + HField tauXRotated, tauYRotated; + rotateVectorFromGreenland(tauX, tauY, tauXRotated, tauYRotated); + + int kinfo; + const int dimension0 = ModelArray::dimensions(ModelArray::Type::H)[0]; + const int dimension1 = ModelArray::dimensions(ModelArray::Type::H)[1]; + + Logged::debug("OASISCoupledOcean: Sending " + TauXKey + " with id " + + std::to_string(couplingId.at(TauXKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(TauXKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &tauXRotated[0], OASIS_No_Restart, &kinfo)); + + Logged::debug("OASISCoupledOcean: Sending " + TauYKey + " with id " + + std::to_string(couplingId.at(TauYKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(TauYKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &tauYRotated[0], OASIS_No_Restart, &kinfo)); + + Logged::debug("OASISCoupledOcean: Sending " + EMPKey + " with id " + + std::to_string(couplingId.at(EMPKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(EMPKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &fwFlux[0], OASIS_No_Restart, &kinfo)); + + const HField qswDown = -qswNet; + Logged::debug("OASISCoupledOcean: Sending " + QSWKey + " with id " + + std::to_string(couplingId.at(QSWKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(QSWKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &qswDown[0], OASIS_No_Restart, &kinfo)); + + const HField qNoSunDown = -qNoSun; + Logged::debug("OASISCoupledOcean: Sending " + QNoSunKey + " with id " + + std::to_string(couplingId.at(QNoSunKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(QNoSunKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &qNoSunDown[0], OASIS_No_Restart, &kinfo)); + + Logged::debug("OASISCoupledOcean: Sending " + SFluxKey + " with id " + + std::to_string(couplingId.at(SFluxKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(SFluxKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &sFlux[0], OASIS_No_Restart, &kinfo)); + + // NEMO wants this field, even if it can be deduced from tauX and tauY + const HField tauMod = (tauX * tauX + tauY * tauY).sqrt(); + Logged::debug("OASISCoupledOcean: Sending " + TauModKey + " with id " + + std::to_string(couplingId.at(TauModKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(TauModKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &tauMod[0], OASIS_No_Restart, &kinfo)); + + // Implicitly copy the 0th DG component (the mean) + const HField cice0 = cice; + Logged::debug("OASISCoupledOcean: Sending " + CIceKey + " with id " + + std::to_string(couplingId.at(CIceKey))); + OASIS_CHECK_ERR(oasis_c_put(couplingId.at(CIceKey), OASISTime, dimension0, dimension1, 1, + OASIS_DOUBLE, OASIS_COL_MAJOR, &cice0[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__) + ": " + NoOASISError); +#endif +} + +void OASISCoupledOcean::configure() +{ + Finalizer::registerUnique(Module::finalize); + Finalizer::registerUnique(Module::finalize); + + 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); + FracQSWKey = Configured::getConfiguration(FracQSWConfigKey, FracQSWKeyDefault); + 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, FracQSWKey }; + if (Configured::getConfiguration(exchangeFirstLayerConfigKey, exchangeFirstLayerDefault)) { + 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(exchangeFirstLayerDefault), "", + "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); +} + +void OASISCoupledOcean::updateTf(const size_t i, const TimestepTime& tst) +{ + tf[i] = Module::getImplementation()(sss[i]); +} +} /* namespace Nextsim */ diff --git a/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp b/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp index 045529ad1..9253c19f0 100644 --- a/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp +++ b/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp @@ -122,6 +122,9 @@ void TOPAZOcean::setData(const ModelState::DataMap& ms) }); slabOcean.setData(ms); + + // We assume that all incoming shortwave radiation is absorbed in the mixed layer + fracQSWAbs = 1.; } void TOPAZOcean::updateTf(size_t i, const TimestepTime& tst) { tf[i] = (*pFreezingPoint)(sss[i]); } diff --git a/physics/src/modules/OceanBoundaryModule/include/OASISCoupledOcean.hpp b/physics/src/modules/OceanBoundaryModule/include/OASISCoupledOcean.hpp new file mode 100644 index 000000000..f3cb97cbf --- /dev/null +++ b/physics/src/modules/OceanBoundaryModule/include/OASISCoupledOcean.hpp @@ -0,0 +1,99 @@ +/*! + * @file OASISCoupledOcean.hpp + * + * @date 21 Mar 2025 + * @author Tim Spain + * @author Einar Ólason + */ + +#ifndef OASISCOUPLEDOCEAN_HPP +#define OASISCOUPLEDOCEAN_HPP + +#include "include/IOceanBoundary.hpp" + +#include "include/Configured.hpp" +#include "include/OASISCoupled.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 exchangeFirstLayerDefault = 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 FracQSWKeyDefault = "I_FrcQsr"; +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_sic"; + +static const std::string SSTConfigKey = moduleName + ".sea_surface_temperature"; +static const std::string SSSConfigKey = moduleName + ".sea_surface_salinity"; +static const std::string UOceanConfigKey = moduleName + ".ocean_u_velocity"; +static const std::string VOceanConfigKey = moduleName + ".ocean_v_velocity"; +static const std::string SSHConfigKey = moduleName + ".sea_surface_height"; +static const std::string FracQSWConfigKey = moduleName + ".frac_solar_absorbed"; +static const std::string MLDConfigKey + = moduleName + ".first_ocean_layer_depth"; // This one is optional +static const std::string TauXConfigKey = moduleName + ".ice_ocean_stress_x"; +static const std::string TauYConfigKey = moduleName + ".ice_ocean_stress_y"; +static const std::string TauModConfigKey = moduleName + ".ice_ocean_stress_modulo"; +static const std::string EMPConfigKey = moduleName + ".fresh_water_flux"; +static const std::string QNoSunConfigKey = moduleName + ".non_solar_heatflux"; +static const std::string QSWConfigKey = moduleName + ".short_wave_flux"; +static const std::string SFluxConfigKey = moduleName + ".salt_flux"; +static const std::string CIceConfigKey = moduleName + ".sea_ice_concentration"; + +//* Ocean boundary data values that are hardcoded. +class OASISCoupledOcean : public IOceanBoundary, + public OASISCoupled, + public Configured { +public: + OASISCoupledOcean() = default; + ~OASISCoupledOcean() override { OASISCoupled::~OASISCoupled(); } + + std::string getName() const override { return moduleName; } + void updateBefore(const TimestepTime& tst) override; + void updateAfter(const TimestepTime& tst) override; + void setData(const ModelState::DataMap& ms) 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; + + void updateTf(size_t i, const TimestepTime& tst); + + // 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, FracQSWKey, MLDKey, TauXKey, TauYKey, + TauModKey, EMPKey, QNoSunKey, QSWKey, SFluxKey, CIceKey; + + // A map relating the namcouple strings to ModelArrayRef's + std::map fieldTags; + + std::vector cplStringsIn; + std::vector cplStringsOut; +}; + +} /* namespace Nextsim */ + +#endif /* OASISCOUPLEDOCEAN_HPP */ diff --git a/physics/src/modules/OceanBoundaryModule/module.cfg b/physics/src/modules/OceanBoundaryModule/module.cfg index ddadaf128..9fa2f8613 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 = true diff --git a/physics/src/modules/include/IAtmosphereBoundary.hpp b/physics/src/modules/include/IAtmosphereBoundary.hpp index 90f3c0d4e..78bfc9f62 100644 --- a/physics/src/modules/include/IAtmosphereBoundary.hpp +++ b/physics/src/modules/include/IAtmosphereBoundary.hpp @@ -5,6 +5,7 @@ #include "include/CheckingModelComponent.hpp" #include "include/ModelArrayRef.hpp" +#include "include/ModelMetadata.hpp" #include "include/Time.hpp" #ifndef IATMOSPHEREBOUNDARY_HPP @@ -52,6 +53,7 @@ class IAtmosphereBoundary : public CheckingModelComponent { virtual ~IAtmosphereBoundary() = default; std::string getName() const override { return "IAtmosphereBoundary"; } + virtual void setMetadata(const ModelMetadata& metadata) { } void setData(const ModelState::DataMap& ms) override { qia.resize(); diff --git a/physics/src/modules/include/IOceanBoundary.hpp b/physics/src/modules/include/IOceanBoundary.hpp index 7aa5bdfb3..67785d88e 100644 --- a/physics/src/modules/include/IOceanBoundary.hpp +++ b/physics/src/modules/include/IOceanBoundary.hpp @@ -7,6 +7,7 @@ #define IOCEANBOUNDARY_HPP #include "include/CheckingModelComponent.hpp" +#include "include/ModelMetadata.hpp" #include "include/constants.hpp" #include "include/gridNames.hpp" @@ -31,6 +32,7 @@ class IOceanBoundary : public CheckingModelComponent { , sFlux(ModelArray::Type::H) , qswow(ModelArray::Type::H, { -1e3, 1e-6 }) , qswBase(ModelArray::Type::H, { -1e3, 1e-6 }) + , fracQSWAbs(ModelArray::Type::H) , tauX(ModelArray::Type::H, { -10, 10 }) , tauY(ModelArray::Type::H, { -10, 10 }) , cice(getStore()) @@ -75,6 +77,7 @@ class IOceanBoundary : public CheckingModelComponent { virtual ~IOceanBoundary() = default; std::string getName() const override { return "IOceanBoundary"; } + virtual void setMetadata(const ModelMetadata& metadata) { } void setData(const ModelState::DataMap& ms) override { qio.resize(); @@ -92,6 +95,7 @@ class IOceanBoundary : public CheckingModelComponent { sFlux.resize(); qswow.resize(); qswBase.resize(); + fracQSWAbs.resize(); tauX.resize(); tauY.resize(); @@ -155,7 +159,7 @@ class IOceanBoundary : public CheckingModelComponent { void mergeFluxesElement(size_t i, const TimestepTime& tst) { // Heat fluxes - partitioned in solar and non-solar - qswNet[i] = cice[i] * qswBase[i] + (1 - cice[i]) * qswow[i]; + qswNet[i] = fracQSWAbs[i] * (cice[i] * qswBase[i] + (1 - cice[i]) * qswow[i]); qNoSun[i] = cice[i] * qio[i] + (1 - cice[i]) * qow[i] - qswNet[i]; // Mass fluxes - fresh water and salt @@ -194,6 +198,7 @@ class IOceanBoundary : public CheckingModelComponent { HField sFlux; // Net surface ocean salt flux, kg m⁻² HField qswow; // Shortwave flux in open water W m⁻² HField qswBase; // Shortwave flux at the base of the ice W m⁻² + HField fracQSWAbs; // The fraction of solar radiation absorbed in the surface ocean layer HField tauX; // x(east)-ward total ocean stress, Pa HField tauY; // y(north)-ward total ocean stress, Pa @@ -201,7 +206,7 @@ class IOceanBoundary : public CheckingModelComponent { ModelArrayRef cice; ModelArrayRef tauXIO; - ModelArrayRef tauYIO; + ModelArrayRef tauYIO; ModelArrayRef evap; ModelArrayRef rain; ModelArrayRef newIce; diff --git a/physics/test/CMakeLists.txt b/physics/test/CMakeLists.txt index af5e6cbe2..cb015578c 100644 --- a/physics/test/CMakeLists.txt +++ b/physics/test/CMakeLists.txt @@ -1,4 +1,19 @@ -if(NOT ENABLE_MPI) +if(ENABLE_MPI) + if (ENABLE_OASIS) + add_executable(testOASISCoupledOcean_MPI1 + "OASISCoupledOcean_test.cpp" "../../core/test/MainMPI.cpp") + target_include_directories(testOASISCoupledOcean_MPI1 + PRIVATE "${ModulesRoot}/OceanBoundaryModule") + target_compile_definitions( + testOASISCoupledOcean_MPI1 + PRIVATE USE_MPI TEST_FILES_DIR=\"${CMAKE_CURRENT_SOURCE_DIR}\" + ) + target_link_libraries(testOASISCoupledOcean_MPI1 PRIVATE nextsimlib doctest::doctest) + + add_test(NAME testOASISCoupledOcean_MPI1 COMMAND + ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 1 testOASISCoupledOcean_MPI1 -ns) + endif () +else() set(COMMON_INCLUDE_DIRS "../../core/src/discontinuousgalerkin" "../../core/src" "../../core/src/modules" "../src" "../src/modules") 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..67afeaded --- /dev/null +++ b/physics/test/OASISCoupledOcean_test.cpp @@ -0,0 +1,125 @@ +/*! + * @file OASISCoupledOcean_test.cpp + * + * @date 15 Feb 2025 + * @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 }); + + double sstIn = -1.84; + double sssIn = 28.0; + // double mldIn = 14.8; + double uIn = -0.14; + double vIn = 0.71; + double sshIn = 14.8; + + HField cice(ModelArray::Type::H); + cice = 0.8; + ModelComponent::getStore().registerArray(Shared::C_ICE_DG, &cice, RO); + + HField tauXIO(ModelArray::Type::H); + tauXIO = -3e-2; + ModelComponent::getStore().registerArray(Protected::IO_STRESS_X, &tauXIO); + + HField tauYIO(ModelArray::Type::H); + tauYIO = 4e-2; + ModelComponent::getStore().registerArray(Protected::IO_STRESS_Y, &tauYIO); + + HField newIce(ModelArray::Type::H); + newIce = 4e-2; + ModelComponent::getStore().registerArray(Shared::NEW_ICE, &newIce, RW); + + HField deltaHice(ModelArray::Type::H); + deltaHice = 1e-4; + ModelComponent::getStore().registerArray(Shared::DELTA_HICE, &deltaHice, RW); + + HField deltaSmelt(ModelArray::Type::H); + deltaSmelt = 1e-4; + ModelComponent::getStore().registerArray(Shared::HSNOW_MELT, &deltaSmelt, RW); + + HField qow(ModelArray::Type::H); + qow = 100; + ModelComponent::getStore().registerArray(Shared::Q_OW, &qow, RW); + + HField tauXOW(ModelArray::Type::H); + tauXOW = tauXIO; + ModelComponent::getStore().registerArray(Shared::OW_STRESS_X, &tauXOW, RW); + + HField tauYOW(ModelArray::Type::H); + tauYOW = tauYIO; + ModelComponent::getStore().registerArray(Shared::OW_STRESS_Y, &tauYOW, RW); + + 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()); + + TimePoint t1("2000-01-01T00:00:00Z"); + TimestepTime tst = { t1, Duration(600) }; + + ocpl.updateBefore(tst); + + ModelArrayRef sst(ModelComponent::getStore()); + ModelArrayRef sss(ModelComponent::getStore()); + ModelArrayRef u(ModelComponent::getStore()); + ModelArrayRef v(ModelComponent::getStore()); + ModelArrayRef ssh(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 ; + std::cout << "Received OCEAN_V at time " << ocpl.OASISTime << ": " << ssh[0] << std::endl ; + */ + REQUIRE(sst[0] == sstIn); + REQUIRE(sss[0] == sssIn); + REQUIRE(u[0] == uIn); + REQUIRE(v[0] == vIn); + REQUIRE(ssh[0] == sshIn); + // REQUIRE(mld[0] == mldIn); + + ModelArrayRef qSwBase(ModelComponent::getStore()); + qSwBase[0] = -2.; + + ModelArrayRef qswow(ModelComponent::getStore()); + qswow[0] = -10.; + + ocpl.updateAfter(tst); + + /* The OASIS output file should contain: + * I_taux: -0.03 + * I_tauy: 0.04 + * I_taumod: 0.05 + * I_fwflux: 0.0609935 + * I_rsso: 3.6 + * I_rsnos: 1651.14 + * I_sfi: 0.000306278 + * I_conc: 0.8 + */ + + 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..54972387e --- /dev/null +++ b/physics/test/generate_OASIS_NetCDF.sh @@ -0,0 +1,75 @@ +#!/bin/bash + +InVars=(I_SST:-1.84 I_SSS:28 I_Uocn:-0.14 I_Vocn:0.71 I_SSH:14.8) +# Optionally add +# InVars+=(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