diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 8c03bc156857..f423372ecd10 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -28,3 +28,69 @@ Omega: Tracers: Base: [Temp, Salt] Debug: [Debug1, Debug2, Debug3] + IOStreams: + # InitialState should only be used when starting from scratch + # After the simulations initial start, the frequency should be + # changed to never so that the initial state file is not read. + InitialState: + UsePointerFile: false + Filename: OmegaMesh.nc + Mode: read + Precision: double + Freq: 1 + FreqUnits: OnStartup + UseStartEnd: false + Contents: + - Restart + # Restarts are used to initialize for all job submissions after the very + # first startup job. We use UseStartEnd with a start time just after the + # simulation start time so that omega does not attempt to use a restart + # for the first startup job. + RestartRead: + UsePointerFile: true + PointerFilename: ocn.pointer + Mode: read + Precision: double + Freq: 1 + FreqUnits: OnStartup + UseStartEnd: true + StartTime: 0001-01-01_00:00:01 + EndTime: 99999-12-31_00:00:00 + Contents: + - Restart + RestartWrite: + UsePointerFile: true + PointerFilename: ocn.pointer + Filename: ocn.restart.$Y-$M-$D_$h.$m.$s + Mode: write + IfExists: replace + Precision: double + Freq: 6 + FreqUnits: months + UseStartEnd: false + Contents: + - Restart + History: + UsePointerFile: false + Filename: ocn.hist.$SimTime + Mode: write + IfExists: replace + Precision: double + Freq: 1 + FreqUnits: months + UseStartEnd: false + Contents: + - Tracers + Highfreq: + UsePointerFile: false + Filename: ocn.hifreq.$Y-$M-$D_$h.$m.$s + Mode: write + IfExists: replace + Precision: single + Freq: 10 + FreqUnits: days + UseStartEnd: true + StartTime: 0001-06-01_00:00:00 + EndTime: 0001-06-30_00:00:00 + Contents: + - Tracers diff --git a/components/omega/doc/devGuide/IOStreams.md b/components/omega/doc/devGuide/IOStreams.md new file mode 100644 index 000000000000..e4174fb6b8e0 --- /dev/null +++ b/components/omega/doc/devGuide/IOStreams.md @@ -0,0 +1,89 @@ +(omega-dev-iostreams)= + +## IO Streams (IOStream) + +Most input and output for Omega occurs through IOStreams. Each stream +defines a file, the contents to be read/written and the time frequency +for reading and writing. Defining streams via the input configuration +file is described in the [User Guide](#omega-user-iostreams). IOStreams +are built on top of the parallel IO infrastructure described in the +[IO Section](#omega-dev-IO) and the Field and Metadata described in the +[Field Section](#omega-dev-Field). Here we describe the classes and functions +used to implement IOStreams. Any module accessing an IOStream instance +or related functions must include the ``IOStream.h`` header file. + +All IOStreams are initialized in a two-step process. A call to the +init routine should take place early in the Omega initialization after +the ModelClock has been initialized using: +```c++ + int Err = IOStream::init(ModelClock); +``` +This routine extracts all the stream definitions from the input configuration +file and creates all the Streams. This initialization also defines the +contents of each Stream but does not yet validate those contents against all +the defined Fields. The contents of all streams should be validated at the +end of initialization (when all Fields have been defined) using the call: +```c++ + bool AllValidate = IOStream::validateAll(); +``` +However, if a stream is needed (eg a read stream) during initialization +before the validateAll call, a single stream can be validated using +```c++ + bool Validated = MyStream.validate(); +``` +and the validation status can be checked with +```c++ + bool Validate = MyStream.isValidated(); +``` +All streams must be validated before use to make sure the Fields have +been defined and the relevant data arrays have been attached to Fields and +are available to access. At the end of a simulation, IOStreams must be +finalized using +```c++ + int Err = IOStream::finalize(ModelClock); +``` +so that any final writes can take place for the OnShutdown streams and to +deallocate all defined streams and arrays. If a stream needs to be removed +before that time, an erase function is provided: +```c++ + IOStream::erase(StreamName); +``` + +For most output streams, we provide a writeAll interface that should be placed +at an appropriate time during the time step loop: +```c++ + int Err = IOStream::writeAll(ModelClock); +``` +This function checks each write stream and writes the file if it is time, based +on a time manager alarm that is defined during initialization for each stream +based on the time frequency in the streams configuration. After writing the +file, the alarm is reset for the next write time. If a file must be written +outside of this routine, a single-stream write can take place using: +```c++ + int Err = IOStream::write(StreamName, ModelClock); +``` + +Reading files (eg for initialization, restart or forcing) does not often +take place all at once, so no readAll interface is provided. Instead, each +input stream is read using: +```c++ + int Err = IOStream::read(StreamName, ModelClock, ReqMetadata); +``` +where ReqMetadata is a variable of type Metadata (defined in Field but +essentially a ``std::map`` for the name/value pair). +This variable should incude the names of global metadata that are desired +from the input file. For example, if a time string is needed to verify the +input file corresponds to a desired time, the required metadata can be +initialized with +```c++ + Metadata ReqMetadata; + ReqMetadata["ForcingTime"] = ""; +``` +The Metadata corresponding to ForcingTime will then be read from the file +and inserted as the Metadata value. If no metadata is to be read from the +file, then an empty ReqMetadata variable can be passed. + +As described in the [User Guide](#omega-user-iostreams), all streams are +defined in the input configuration file and most other IOStream functions +are associated either with that initialization or to support the read/write +functions above. diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 240c51e70f69..711853fc6f6f 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -32,6 +32,7 @@ userGuide/Decomp userGuide/Dimension userGuide/Field userGuide/IO +userGuide/IOStreams userGuide/Halo userGuide/HorzMesh userGuide/HorzOperators @@ -66,6 +67,7 @@ devGuide/Decomp devGuide/Dimension devGuide/Field devGuide/IO +devGuide/IOStreams devGuide/Halo devGuide/HorzMesh devGuide/HorzOperators diff --git a/components/omega/doc/userGuide/IOStreams.md b/components/omega/doc/userGuide/IOStreams.md new file mode 100644 index 000000000000..475a0fc0d3ca --- /dev/null +++ b/components/omega/doc/userGuide/IOStreams.md @@ -0,0 +1,142 @@ +(omega-user-iostreams)= + +## IO Streams (IOStream) + +IO Streams are the primary mechanism for users to specify input and output +for Omega. An IOStream can be defined for any number of fields and at desired +time frequencies (including one-time or at startup/shutdown). IOStreams are +defined in the Omega input configuration file in an IOStreams section: + +```yaml +Omega: + # other config options removed for brevity + IOStreams: + InitialState: + UsePointerFile: false + Filename: OmegaMesh.nc + Mode: read + Precision: double + Freq: 1 + FreqUnits: OnStartup + UseStartEnd: false + Contents: + - Restart + RestartWrite: + UsePointerFile: true + PointerFilename: ocn.pointer + Filename: ocn.restart.$Y-$M-$D_$h.$m.$s + Mode: write + IfExists: replace + Precision: double + Freq: 6 + FreqUnits: months + UseStartEnd: false + Contents: + - Restart + History: + UsePointerFile: false + Filename: ocn.hist.$SimTime + Mode: write + IfExists: replace + Precision: double + Freq: 1 + FreqUnits: months + UseStartEnd: false + Contents: + - Tracers + Highfreq: + UsePointerFile: false + Filename: ocn.hifreq.$Y-$M-$D_$h.$m.$s + Mode: write + IfExists: replace + Precision: single + Freq: 10 + FreqUnits: days + UseStartEnd: true + StartTime: 0001-06-01_00.00.00 + EndTime: 0001-06-30_00.00.00 + Contents: + - Tracers +``` + +Each stream has a number of required and optional parameters for customizing +input and output. These options are indented below the stream name as shown +in the sample YAML entries above. They include: +- **UsePointerFile:** A required flag that is either true or false. A pointer +file is used for cases like restart files where the last file written can +be stored for the next job submission so that the configuration file does +not need to be edited between job submissions. +- **PointerFilename:** Only required if UsePointerFile is true and should +be set to the full filename (with path) for the pointer file. Each stream +using a pointer file must define a unique pointer file name. +- **Filename:** Required in all cases except input streams using a pointer +file. This is the complete name (with path) of the file to be read or written. +A filename template is also supported in which simulation (or real) time +can be used in the file name. As the examples above show, accepted keys for +a template can be: + - $SimTime for the current simulation time in a standard time string (note + that this time string may include colon separators that can be a problem + for filenames so using the individual keys below is preferred). + - $Y for the current simulation year + - $M for the current simulation month + - $D for the current simulation day + - $h for the current simulation hour + - $m for the current simulation minute + - $s for the current simulation second + - $WallTime for the time IRL for use when you might need the actual time for + a debug time stamp +- **Mode:** A required field that is either read or write. There is no + readwrite option (eg for restarts) so a separate stream should be + defined for such cases as in the examples above. +- **IfExists:** A required field for write streams that determines behavior + if the file already exists. Acceptable options are: + - Fail if you want the code to exit with an error + - Replace if you want to replace the existing file with the new file + - Append if you want to append (eg multiple time slices) to the existing + file (this option is not currently supported). +- **Precision:** A field that determines whether floating point numbers are + written in full (double) precision or reduced (single). Acceptable values + are double or single. If not present, double is assumed, but a warning + message will be generated so it is best to explicitly include it. +- **Freq:** A required integer field that determines the frequency of + input/output in units determined by the next FreqUnits entry. +- **FreqUnits:** A required field that, combined with the integer frequency, + determines the frequency of input/output. Acceptable values include: + - OnStartup for files read/written once on startup + - OnShutdown for files read/written only once on model exit + - AtTime or OnTime or Time or TimeInstant for a one-time read or write + at the time specified in the StartTime entry + - Never if the stream should not be used but you wish to retain the + entry in the config file (a warning will be output to log file) + - Years for a frequency every Freq years (*not* Freq times per year) + - Months for a frequency every Freq months (*not* Freq times per month) + - Days for a frequency every Freq days (*not* Freq times per day) + - Hours for a frequency every Freq hours (*not* Freq times per hour) + - Minutes for a frequency every Freq minutes (*not* Freq times per minute) + - Seconds for a frequency every Freq seconds (*not* Freq times per seconds) +- **UseStartEnd:** A required entry that is true or false and is used if the + I/O is desired only within a certain time interval. An example might be + for specifying high-frequency output within a certain period of a simulation. +- **StartTime:** A field only required when UseStartEnd is true or if + the FreqUnits request a one-time read/write. The StartTime must be a time + string of the format YYYY-MM-DD_hh.mm.ss (though the delimiters can be + any non-numeric character). The year entry is the integer year and can be + four or more digits. The StartTime is inclusive - the I/O will occur at or + after that date/time. +- **EndTime:** A field that is only required when UseStartEnd is true. It + requires the same format as StartTime but unlike StartTime, the EndTime + is not inclusive and I/O only occurs for times before the EndTime. If a + file is desired at the EndTime, the user should specify an EndTime slightly + later (less than a time step) than the desired end time. +- **Contents:** This is a required field that contains an itemized list of + each Field or FieldGroup that is desired in the output. The name must + match a name of a defined Field or Group within Omega. Group names are + preferred to keep the list of fields short so Omega will define convenient + FieldGroups like Restart, State, Tracers that will include all members + of the group. If only a subset of Fields from a Group is desired, the + individual Field names should be specified and not the Group name. + +This streams configuration should be sufficient to define all input and output +from the model and provides a relatively simple interface for a typical user. +However, if necessary (eg before streams have been defined), the specific +interfaces in the lower level [IO](#omega-user-IO) module can be used. diff --git a/components/omega/src/base/IO.cpp b/components/omega/src/base/IO.cpp index 25aa988a8ed4..c47d7564e279 100644 --- a/components/omega/src/base/IO.cpp +++ b/components/omega/src/base/IO.cpp @@ -297,8 +297,18 @@ int openFile( // Closes an open file using the fileID, returns an error code int closeFile(int &FileID /// [in] ID of the file to be closed ) { - // Just calls the PIO close routine - int Err = PIOc_closefile(FileID); + int Err = 0; + + // Make sure all operations completed before closing + Err = PIOc_sync(FileID); + if (Err != PIO_NOERR) + LOG_ERROR("Error syncing file before closing"); + + // Call the PIO close routine + Err = PIOc_closefile(FileID); + if (Err != PIO_NOERR) + LOG_ERROR("Error closing file {} in PIO", FileID); + return Err; } // End closeFile @@ -454,9 +464,9 @@ int writeMeta(const std::string &MetaName, // [in] name of metadata } // End writeMeta (R8) -int writeMeta(const std::string &MetaName, // [in] name of metadata - std::string MetaValue, // [in] value of metadata - int FileID, // [in] ID of the file for writing +int writeMeta(const std::string &MetaName, // [in] name of metadata + const std::string &MetaValue, // [in] value of metadata + int FileID, // [in] ID of the file for writing int VarID // [in] ID for variable associated with metadata ) { @@ -690,8 +700,10 @@ int readArray(void *Array, // [out] array to be read // Find variable ID from file Err = PIOc_inq_varid(FileID, VarName.c_str(), &VarID); - if (Err != PIO_NOERR) + if (Err != PIO_NOERR) { LOG_ERROR("IO::readArray: Error finding varid for variable {}", VarName); + return Err; + } // PIO Read array call to read the distributed array PIO_Offset ASize = Size; @@ -721,6 +733,19 @@ int writeArray(void *Array, // [in] array to be written PIO_Offset Asize = Size; Err = PIOc_write_darray(FileID, VarID, DecompID, Asize, Array, FillValue); + if (Err != PIO_NOERR) { + LOG_ERROR("Error in PIO writing distributed array"); + return Err; + } + + // Make sure write is complete before returning + // We may be able to remove this for efficiency later but it was + // needed during testing + Err = PIOc_sync(FileID); + if (Err != PIO_NOERR) { + LOG_ERROR("Error in PIO sychronizing file after write"); + return Err; + } return Err; diff --git a/components/omega/src/base/IO.h b/components/omega/src/base/IO.h index 7c61c47c08dd..8472fee92a65 100644 --- a/components/omega/src/base/IO.h +++ b/components/omega/src/base/IO.h @@ -194,9 +194,9 @@ int writeMeta(const std::string &MetaName, ///< [in] name of metadata int FileID, ///< [in] ID of the file for writing int VarID ///< [in] ID for variable associated with metadata ); -int writeMeta(const std::string &MetaName, ///< [in] name of metadata - std::string MetaValue, ///< [in] value of metadata - int FileID, ///< [in] ID of the file for writing +int writeMeta(const std::string &MetaName, ///< [in] name of metadata + const std::string &MetaValue, ///< [in] value of metadata + int FileID, ///< [in] ID of the file for writing int VarID ///< [in] ID for variable associated with metadata ); diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp new file mode 100644 index 000000000000..ebbf0c2d9105 --- /dev/null +++ b/components/omega/src/infra/IOStream.cpp @@ -0,0 +1,2799 @@ +//===-- infra/IOStream.cpp - IO stream implementation -----------*- C++ -*-===// +// +// This file implements classes and methods for IO Streams. IOStreams define +// reading and writing of fields from/to a data file. Each stream contains +// information on the file location, the frequency of input/output, the +// fields to be read/written and other information necessary for file I/O. +// Note that this is different from the C++ stdlib iostreams +// +//===----------------------------------------------------------------------===// + +#include "IOStream.h" +#include "DataTypes.h" +#include "Dimension.h" +#include "Field.h" +#include "IO.h" +#include "Logging.h" +#include "OmegaKokkos.h" +#include "TimeMgr.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace OMEGA { + +// Create static class members +std::map> IOStream::AllStreams; + +//------------------------------------------------------------------------------ +// Initializes all streams defined in the input configuration file. This +// does not validate the contents of the streams since the relevant Fields +// may not have been defined yet. Returns an error code. +int IOStream::init(Clock &ModelClock //< [inout] Omega model clock +) { + + int Err = 0; // default return code + + // Retrieve the model configuration and get the streams sub-config + Config *OmegaConfig = Config::getOmegaConfig(); + Config StreamsCfgAll("IOStreams"); + Err = OmegaConfig->get(StreamsCfgAll); + if (Err != 0) { // could not find the streams configuration + LOG_ERROR("Could not find Streams configuration in Omega input file"); + return Err; + } + + // Loop over all streams in the subconfiguration and create them + for (auto It = StreamsCfgAll.begin(); It != StreamsCfgAll.end(); ++It) { + + // Get the stream name and the sub-configuration associated with it + std::string StreamName = It->first.as(); + Config StreamCfg(StreamName); + Err = StreamsCfgAll.get(StreamCfg); + if (Err != 0) { + LOG_ERROR("Error retrieving configuration for stream {}", StreamName); + return Err; + } + + // Call the create routine to create the stream + Err = create(StreamName, StreamCfg, ModelClock); + if (Err != 0) { + LOG_ERROR("Error creating stream {}", StreamName); + return Err; + } + + } // end loop over all streams + + return Err; + +} // End initialize + +//------------------------------------------------------------------------------ +// Performs a final write of any streams that have the OnShutdown option and +// then removes all streams to clean up. Returns an error code. +int IOStream::finalize( + const Clock &ModelClock // [in] Model clock needed for time stamps +) { + + int Err = 0; + bool FinalCall = true; + + // Loop over all streams and call write function for any write streams + // with the OnShutdown flag + for (auto Iter = AllStreams.begin(); Iter != AllStreams.end(); Iter++) { + + std::string StreamName = Iter->first; + std::shared_ptr ThisStream = Iter->second; + bool ForceWrite = false; + + int Err1 = 0; + if (ThisStream->OnShutdown) + Err1 = ThisStream->writeStream(ModelClock, ForceWrite, FinalCall); + + if (Err1 != 0) { + LOG_ERROR("Error trying to write stream {} at shutdown", StreamName); + ++Err; + } + } // end loop over streams + + // Remove all streams + AllStreams.clear(); + + return Err; + +} // End finalize + +//------------------------------------------------------------------------------ +// Retrieves a pointer to a previously defined stream. +std::shared_ptr +IOStream::get(const std::string &StreamName ///< [in] name of stream to retrieve +) { + // Find stream in list of streams and return the pointer + if (AllStreams.find(StreamName) != AllStreams.end()) { + return AllStreams[StreamName]; + } else { + LOG_ERROR("Cannot retrieve Stream {}. Stream has not been created.", + StreamName); + return nullptr; + } + +} // End get stream + +//------------------------------------------------------------------------------ +// Adds a field to the contents of a stream. Because streams may be created +// before all Fields have been defined, we only store the name. Validity +// is either checked during read/write or can be checked using the validate +// function. +void IOStream::addField(const std::string &FieldName ///< [in] Name of field +) { + this->Contents.insert(FieldName); +} // End addField + +//------------------------------------------------------------------------------ +// Removes a field from the contents. Provided for symmetry, but not +// typically used. +void IOStream::removeField(const std::string &FieldName ///< [in] Name of field +) { + this->Contents.erase(FieldName); +} // End removeField + +//------------------------------------------------------------------------------ +// Validates an IOStream and its contents. This must be called after all +// relevant fields have been defined, typically toward the end of omega model +// initialization. The routine also expands all group names to the individual +// field names that are members of the group. Returns true if all contents and +// variables are valid. +bool IOStream::validate() { + + bool ReturnVal = true; + // Return if already validated + if (Validated) + return ReturnVal; + + // Expand group names to list of individual fields + // First identify any group names in the Contents + std::set GroupNames; + for (auto IField = Contents.begin(); IField != Contents.end(); ++IField) { + std::string FieldName = *IField; + // If this name is a group name, add it to the list + if (FieldGroup::exists(FieldName)) + GroupNames.insert(FieldName); + } + + // Now for each group, extract field names and add to contents if it + // is not there already. Since Contents is a set we can just insert + // and it will take care of duplicate fields. + for (auto IGrp = GroupNames.begin(); IGrp != GroupNames.end(); ++IGrp) { + std::string GrpName = *IGrp; + // Get Fields from Group + std::set FieldList = + FieldGroup::getFieldListFromGroup(GrpName); + // Insert field names into Contents + for (auto IField = FieldList.begin(); IField != FieldList.end(); + ++IField) { + Contents.insert(*IField); + } + // Remove group name from contents + Contents.erase(GrpName); + } + + // Loop through all the field names in Contents and check whether they + // have been defined as an Field + for (auto IField = Contents.begin(); IField != Contents.end(); ++IField) { + std::string FieldName = *IField; + + if (not Field::exists(FieldName)) { + LOG_ERROR("Cannot validate stream {}: Field {} has not been defined", + Name, FieldName); + ReturnVal = false; + } + } + + if (ReturnVal) + Validated = true; + return ReturnVal; + +} // End validate + +//------------------------------------------------------------------------------ +// Checks that a stream has been validated +bool IOStream::isValidated() { return Validated; } + +//------------------------------------------------------------------------------ +// Validates all streams and their contents. If used, this must be called at +// the end of initialization to ensure all Fields have been defined. +// Returns true if all streams are valid. +bool IOStream::validateAll() { + + bool ReturnVal = true; // default is all valid + + // Loop over all streams and call validate function + for (auto Iter = AllStreams.begin(); Iter != AllStreams.end(); Iter++) { + std::string StreamName = Iter->first; + std::shared_ptr ThisStream = Iter->second; + bool Valid = ThisStream->validate(); + if (!Valid) { + ReturnVal = false; + LOG_ERROR("IOStream validateAll: stream {} has invalid entries", + StreamName); + } + } + + return ReturnVal; + +} // End validateAll + +//------------------------------------------------------------------------------ +// Reads a single stream if it is time. Returns an error code. +int IOStream::read( + const std::string &StreamName, // [in] Name of stream + const Clock &ModelClock, // [in] Model clock for time info + Metadata &ReqMetadata, // [inout] global metadata requested from file + bool ForceRead // [in] optional: read even if not time +) { + int Err = 0; // default return code + + // Retrieve stream by name and make sure it exists + auto StreamItr = AllStreams.find(StreamName); + if (StreamItr != AllStreams.end()) { + // Stream found, call the read function + std::shared_ptr ThisStream = StreamItr->second; + Err = ThisStream->readStream(ModelClock, ReqMetadata, ForceRead); + } else { // Stream not found, return error + LOG_ERROR("Unable to read stream {}. Stream not defined", StreamName); + Err = 1; + } + + return Err; + +} // End read stream + +//------------------------------------------------------------------------------ +// Writes a single stream if it is time. Returns an error code. +int IOStream::write( + const std::string &StreamName, // [in] Name of stream + const Clock &ModelClock, // [in] Model clock needed for time stamps + bool ForceWrite // [in] optional: write even if not time +) { + int Err = 0; // default return code + + // Retrieve stream by name and make sure it exists + auto StreamItr = AllStreams.find(StreamName); + if (StreamItr != AllStreams.end()) { + // Stream found, call the write function + std::shared_ptr ThisStream = StreamItr->second; + Err = ThisStream->writeStream(ModelClock, ForceWrite); + } else { + // Stream not found, return error + LOG_ERROR("Unable to write stream {}. Stream not defined", StreamName); + Err = 1; + } + + return Err; + +} // End write stream + +//------------------------------------------------------------------------------ +// Loops through all streams and writes them if it is time. This is +// useful if most I/O is consolidated at one point (eg end of step). +int IOStream::writeAll( + const Clock &ModelClock // [in] Model clock needed for time stamps +) { + + int Err = 0; // accumulated error for return value + + // Loop over all streams and call write function for any write streams + for (auto Iter = AllStreams.begin(); Iter != AllStreams.end(); Iter++) { + + std::string StreamName = Iter->first; + std::shared_ptr ThisStream = Iter->second; + + int Err1; + if (ThisStream->Mode == IO::ModeWrite) { + Err1 = ThisStream->writeStream(ModelClock); + } + + // check for errors + if (Err1 != 0) { + ++Err; + LOG_ERROR("writeAll error in stream {}", StreamName); + } + } + + return Err; + +} // End writeAll + +//------------------------------------------------------------------------------ +// Constructs an empty IOStream +IOStream::IOStream() { + + Name = "Unknown"; + Filename = "Unknown"; + FilenameIsTemplate = false; + ExistAction = IO::IfExists::Fail; + Mode = IO::Mode::ModeUnknown; + ReducePrecision = false; + OnStartup = false; + OnShutdown = false; + UsePointer = false; + PtrFilename = " "; + UseStartEnd = false; + Validated = false; +} + +//------------------------------------------------------------------------------ +// Creates a new stream and adds to the list of all streams, based on +// options in the input model configuration. This routine is called by +// the IOStreams initialize function. It requires an initialized model +// clock and stream alarms are attached to this clock during creation. + +int IOStream::create(const std::string &StreamName, //< [in] name of stream + Config &StreamConfig, //< [in] input stream configuration + Clock &ModelClock //< [inout] Omega model clock +) { + + int Err = 0; + + // Check whether the stream already exists + if (AllStreams.find(StreamName) != AllStreams.end()) { + // Stream already exists, return error + Err = 1; + LOG_ERROR("Attempt to create stream {} that already exists", StreamName); + return Err; + } + + // Create a new pointer and set name + auto NewStream = std::make_shared(); + NewStream->Name = StreamName; + + // Set file mode (Read/Write) + std::string StreamMode; + Err = StreamConfig.get("Mode", StreamMode); + NewStream->Mode = IO::ModeFromString(StreamMode); + if ((Err != 0) or (NewStream->Mode == IO::ModeUnknown)) { + LOG_ERROR("Bad or non-existent Mode for stream {}", StreamName); + Err = 2; + return Err; + } + + // Set file name + // First check if using a pointer file + Err = StreamConfig.get("UsePointerFile", NewStream->UsePointer); + if (Err != 0) { + LOG_ERROR("UsePointerFile flag not found in stream {}", StreamName); + return Err; + } + // If the stream uses a pointer file, get the pointer filename + if (NewStream->UsePointer) { + NewStream->PtrFilename = ""; // initialize to blank filename + Err = StreamConfig.get("PointerFilename", NewStream->PtrFilename); + if (Err != 0) { // pointer filename not found + LOG_ERROR("Pointer filename not found for stream {}", StreamName); + return Err; + } + } + + // For file reads that are using a pointer file, the filename is + // read later from that pointer file. All other cases need to read + // a filename or filename template from config + // Otherwise, read filename from config + std::string StreamFilename = ""; + if (NewStream->Mode != IO::ModeRead or !(NewStream->UsePointer)) { + Err = StreamConfig.get("Filename", StreamFilename); + if (Err != 0) { + LOG_ERROR("Error getting Filename for stream {} from Config", + StreamName); + return Err; + } + // Check to see if filename is a template and needs to be constructed + // later with time information + auto TemplateFound = StreamFilename.find("$"); + if (TemplateFound != std::string::npos) { + NewStream->FilenameIsTemplate = true; + } else { + NewStream->FilenameIsTemplate = false; + } + } + // Add filename to stream + NewStream->Filename = StreamFilename; + + // Set flag to reduce precision for double precision reals. If no flag + // present, assume full (double) precision + std::string PrecisionString; + Err = StreamConfig.get("Precision", PrecisionString); + if (Err != 0) + PrecisionString = "double"; + NewStream->setPrecisionFlag(PrecisionString); + + // Set the action to take if a file already exists + // This is only needed for writes so only perform check for write mode + NewStream->ExistAction = IO::IfExists::Fail; // default is to fail + if (NewStream->Mode == IO::ModeWrite) { + std::string ExistAct; + Err = StreamConfig.get("IfExists", ExistAct); + if (Err == 0) + NewStream->ExistAction = IO::IfExistsFromString(ExistAct); + } + + // Set alarm based on read/write frequency + // Use stream name as alarm name + std::string AlarmName = StreamName; + + // For alarms, need to retrieve clock start time + TimeInstant ClockStart = ModelClock.getStartTime(); + Calendar *CalendarPtr; + Err = ClockStart.get(CalendarPtr); + if (Err != 0) { + LOG_ERROR("Unable to retrieve clock info while constructing stream {}", + StreamName); + return Err; + } + + // Read frequency of input/output + int IOFreq; + Err = StreamConfig.get("Freq", IOFreq); + if (Err != 0) { + LOG_ERROR("Frequency missing for stream {}", StreamName); + return Err; + } + if (IOFreq < 1) { + LOG_ERROR("Invalid frequency {} for IO stream {} ", StreamName); + Err = 3; + return Err; + } + std::string IOFreqUnits; + Err = StreamConfig.get("FreqUnits", IOFreqUnits); + if (Err != 0) { + LOG_ERROR("FreqUnits missing for stream {}", StreamName); + return Err; + } + + // convert string to lower case for easier comparison + std::transform(IOFreqUnits.begin(), IOFreqUnits.end(), IOFreqUnits.begin(), + [](unsigned char C) { return std::tolower(C); }); + + // Based on input frequency and units, create the alarm or set flags + bool HasAlarm = false; + NewStream->OnStartup = false; + NewStream->OnShutdown = false; + + if (IOFreqUnits == "years") { + + TimeInterval AlarmInt(IOFreq, TimeUnits::Years); + NewStream->MyAlarm = Alarm(AlarmName, AlarmInt, ClockStart); + HasAlarm = true; + + } else if (IOFreqUnits == "months") { + + TimeInterval AlarmInt(IOFreq, TimeUnits::Months); + NewStream->MyAlarm = Alarm(AlarmName, AlarmInt, ClockStart); + HasAlarm = true; + + } else if (IOFreqUnits == "days") { + + TimeInterval AlarmInt(IOFreq, TimeUnits::Days); + NewStream->MyAlarm = Alarm(AlarmName, AlarmInt, ClockStart); + HasAlarm = true; + + } else if (IOFreqUnits == "hours") { + + TimeInterval AlarmInt(IOFreq, TimeUnits::Hours); + NewStream->MyAlarm = Alarm(AlarmName, AlarmInt, ClockStart); + HasAlarm = true; + + } else if (IOFreqUnits == "minutes") { + + TimeInterval AlarmInt(IOFreq, TimeUnits::Minutes); + NewStream->MyAlarm = Alarm(AlarmName, AlarmInt, ClockStart); + HasAlarm = true; + + } else if (IOFreqUnits == "seconds") { + + TimeInterval AlarmInt(IOFreq, TimeUnits::Seconds); + NewStream->MyAlarm = Alarm(AlarmName, AlarmInt, ClockStart); + HasAlarm = true; + + } else if (IOFreqUnits == "onstartup") { + + NewStream->OnStartup = true; + + } else if (IOFreqUnits == "onshutdown") { + NewStream->OnShutdown = true; + + } else if (IOFreqUnits == "attime" or IOFreqUnits == "ontime" or + IOFreqUnits == "time" or IOFreqUnits == "timeinstant") { + + // A one-time event for this stream - use the StartTime string + // as the time instant to use + std::string StrtTime; + Err = StreamConfig.get("StartTime", StrtTime); + if (Err == 0) { + TimeInstant AlarmTime(CalendarPtr, StrtTime); + NewStream->MyAlarm = Alarm(AlarmName, AlarmTime); + HasAlarm = true; + } else { + LOG_ERROR("Stream {} requests a one-time read/write but StartTime" + "not provided", + StreamName); + return Err; + } + + } else if (IOFreqUnits == "never") { + + LOG_WARN("Stream {} has IO frequency of never and will be skipped", + StreamName); + Err = 0; + return Err; + + } else { + + if (!NewStream->OnStartup and !NewStream->OnShutdown) { + LOG_ERROR("Unknown IOFreqUnits option for stream {}", StreamName); + Err = 4; + return Err; + } + } + // If an alarm is set, attach it to the model clock + if (HasAlarm) { + Err = ModelClock.attachAlarm(&(NewStream->MyAlarm)); + if (Err != 0) { + LOG_ERROR("Error attaching alarm to model clock for stream {}", + StreamName); + return Err; + } + } + + // Use a start and end time to define an interval in which stream is active + Err = StreamConfig.get("UseStartEnd", NewStream->UseStartEnd); + if (Err != 0) { // Start end flag not in config, assume false + NewStream->UseStartEnd = false; + Err = 0; + } + + // Set Alarms for start and end time + if (NewStream->UseStartEnd) { + std::string StartTimeStr; + std::string EndTimeStr; + int Err = StreamConfig.get("StartTime", StartTimeStr); + if (Err != 0) { + LOG_ERROR("Stream {} requests UseStartEnd but no start time provided", + StreamName); + return Err; + } + Err = StreamConfig.get("EndTime", EndTimeStr); + if (Err != 0) { + LOG_ERROR("Stream {} requests UseStartEnd but no end time provided", + StreamName); + return Err; + } + TimeInstant Start(CalendarPtr, StartTimeStr); + TimeInstant End(CalendarPtr, EndTimeStr); + std::string StartName = StreamName + "Start"; + std::string EndName = StreamName + "End"; + NewStream->StartAlarm = Alarm(StartName, Start); + NewStream->EndAlarm = Alarm(EndName, End); + Err = ModelClock.attachAlarm(&(NewStream->StartAlarm)); + if (Err != 0) { + LOG_ERROR("Error attaching start alarm to model clock for stream {}", + StreamName); + return Err; + } + Err = ModelClock.attachAlarm(&(NewStream->EndAlarm)); + if (Err != 0) { + LOG_ERROR("Error attaching end alarm to model clock for stream {}", + StreamName); + return Err; + } + } // endif UseStartEnd + + // Now we add the list of field names to the stream + // First get the contents list + std::vector FieldContents; + Err = StreamConfig.get("Contents", FieldContents); + if (Err != 0) { + LOG_ERROR("Can not find contents for stream {}", StreamName); + return Err; + } + + // The contents are stored as an ordered set so we use the addField + // interface to add each name. Note that in this context, the field + // name can also be a group name. Group names are expanded during the + // validate stage. + for (int IField = 0; IField < FieldContents.size(); ++IField) { + NewStream->addField(FieldContents[IField]); + } + + // The contents list has not yet been validated. + NewStream->Validated = false; + + // If we have made it to this point, we have a valid stream to add to + // the list + AllStreams[StreamName] = NewStream; + + return Err; + +} // End IOStream create + +//------------------------------------------------------------------------------ +// Define all dimensions used. Returns an error code as well as a map +// of dimension names to defined dimension IDs. +int IOStream::defineAllDims( + int FileID, ///< [in] id assigned to the IO file + std::map &AllDimIDs ///< [out] dim name, assigned ID +) { + + int Err = 0; + + for (auto IDim = Dimension::begin(); IDim != Dimension::end(); ++IDim) { + std::string DimName = IDim->first; + // For back compatibility, we also allow an older name (typically + // with first name lower case). MaxCellsOnEdge is an exception since + // it was named TWO in MPAS + std::string OldDimName = DimName; + if (DimName == "MaxCellsOnEdge") { + OldDimName = "TWO"; + } else { + OldDimName[0] = std::tolower(OldDimName[0]); + } + I4 Length = IDim->second->getLengthGlobal(); + I4 DimID; + + // For input files, we read the DimID from the file + if (Mode == IO::ModeRead) { + // If dimension not found, only generate a warning since there + // may be some dimensions that are not required + I4 InLength; + Err = IO::getDimFromFile(FileID, DimName, DimID, InLength); + if (Err != 0) { // can't find dim in file + // Try again using old name for back compatibility to MPAS + Err = IO::getDimFromFile(FileID, OldDimName, DimID, InLength); + if (Err == 0) { + LOG_INFO("Ignore PIO Error for Dimension {}: ", DimName); + LOG_INFO("Found under old dimension name {}: ", OldDimName); + } else { + if (Err != 0) + LOG_WARN("Dimension {} not found in input stream {}", DimName, + Name); + } + continue; + } + // Check dimension length in input file matches what is expected + if (InLength != Length) { + LOG_ERROR("Inconsistent length for dimension {} in input stream {}", + DimName, Name); + Err = 3; + return Err; + } + } // end read case + + // For output files, we need to define the dimension + if (Mode == IO::ModeWrite) { + + Err = IO::defineDim(FileID, DimName, Length, DimID); + if (Err != 0) { + LOG_ERROR("Error defining dimension {} for output stream {}", + DimName, Name); + return Err; + } + } // end write case + + // Add the DimID to map for later use + AllDimIDs[DimName] = DimID; + + } // end loop over all dims + + return Err; + +} // End defineAllDims + +//------------------------------------------------------------------------------ +// Computes the parallel decomposition (offsets) for a field needed for parallel +// I/O. Return error code and also Decomp ID and array size for field. +int IOStream::computeDecomp( + std::shared_ptr FieldPtr, // [in] pointer to Field + std::map &AllDimIDs, // [in] dimension IDs + int &DecompID, // [out] ID assigned to the decomposition + int &LocalSize, // [out] size of local array + std::vector &DimLengths // [out] vector of local dim lengths +) { + + int Err = 0; + + // Retrieve some basic field information + std::string FieldName = FieldPtr->getName(); + IO::IODataType MyIOType = getFieldIOType(FieldPtr); + int NDims = FieldPtr->getNumDims(); + if (NDims < 1) { + LOG_ERROR("Invalid number of dimensions for Field {}", FieldName); + Err = 1; + return Err; + } + + std::vector DimNames(NDims); + Err = FieldPtr->getDimNames(DimNames); + if (Err != 0) { + LOG_ERROR("Error retrieving dimension names for Field {}", FieldName); + return Err; + } + + // Get dimension and size information for each dimension + // For the offset calculation, we also create dimension and dimension offsets + // in index order up to the max of 5 dimensions. The active indices are + // at the end in array index order. + I4 GlobalSize = 1; + LocalSize = 1; + constexpr I4 MaxDims = 5; + std::vector DimLengthsGlob(NDims, 1); + std::vector DimOffsets(MaxDims); + std::vector DimLengthGlob(MaxDims, 1); // lengths padded to MaxDims + std::vector DimLengthLoc(MaxDims, 1); // lengths padded to MaxDims + + for (int IDim = 0; IDim < NDims; ++IDim) { + I4 StartDim = MaxDims - NDims; + std::string DimName = DimNames[IDim]; + std::shared_ptr ThisDim = Dimension::get(DimName); + DimLengths[IDim] = ThisDim->getLengthLocal(); + DimLengthsGlob[IDim] = ThisDim->getLengthGlobal(); + DimOffsets[StartDim + IDim] = ThisDim->getOffset(); + DimLengthLoc[StartDim + IDim] = DimLengths[IDim]; + DimLengthGlob[StartDim + IDim] = DimLengthsGlob[IDim]; + LocalSize *= DimLengths[IDim]; + GlobalSize *= DimLengthsGlob[IDim]; + } + + // Create the data decomposition based on dimension information + // Compute offset index (0-based global index of each element) in + // linear address space. Needed for the decomposition definition. + // -1 is used to denote indices that are not written + // The linear offset along each dimension has already been computed + // by the dimension class. + + std::vector Offset(LocalSize, -1); + + // Compute strides in linear space for each dimension + // For the proper stride in storage order, we start from the rightmost + // index in index order. + std::vector StrideGlob(MaxDims, 1); + std::vector StrideLoc(MaxDims, 1); + for (int IDim = MaxDims - 2; IDim >= 0; --IDim) { + StrideGlob[IDim] = StrideGlob[IDim + 1] * DimLengthGlob[IDim + 1]; + StrideLoc[IDim] = StrideLoc[IDim + 1] * DimLengthLoc[IDim + 1]; + } + + // Compute full array offsets based on each dimensions linear offset + // Skip -1 entries that are outside owned domain + for (int N = 0; N < DimLengthLoc[0]; ++N) { + I4 NGlob = 0; + if (NDims >= 5) + NGlob = (DimOffsets[0])(N); + if (NGlob < 0) + continue; // index outside owned domain + for (int M = 0; M < DimLengthLoc[1]; ++M) { + I4 MGlob = 0; + if (NDims >= 4) + MGlob = (DimOffsets[1])(M); + if (MGlob < 0) + continue; + for (int K = 0; K < DimLengthLoc[2]; ++K) { + I4 KGlob = 0; + if (NDims >= 3) + KGlob = (DimOffsets[2])(K); + if (KGlob < 0) + continue; + for (int J = 0; J < DimLengthLoc[3]; ++J) { + I4 JGlob = 0; + if (NDims >= 2) + JGlob = (DimOffsets[3])(J); + if (JGlob < 0) + continue; + for (int I = 0; I < DimLengthLoc[4]; ++I) { + I4 IGlob = (DimOffsets[4])(I); + if (IGlob < 0) + continue; + int Add = I * StrideLoc[4] + J * StrideLoc[3] + + K * StrideLoc[2] + M * StrideLoc[1] + + N * StrideLoc[0]; + Offset[Add] = IGlob * StrideGlob[4] + JGlob * StrideGlob[3] + + KGlob * StrideGlob[2] + MGlob * StrideGlob[1] + + NGlob * StrideGlob[0]; + } + } + } + } + } + + Err = OMEGA::IO::createDecomp(DecompID, MyIOType, NDims, DimLengthsGlob, + LocalSize, Offset, OMEGA::IO::DefaultRearr); + if (Err != 0) { + LOG_ERROR("Error creating decomp for field {} in stream {}", FieldName, + Name); + return Err; + } + + return Err; + +} // End computeDecomp + +//------------------------------------------------------------------------------ +// Write all metadata associated with a field +int IOStream::writeFieldMeta( + const std::string FieldName, // [in] metadata for field + int FileID, // [in] id assigned to open file + int FieldID // [in] id assigned to the field +) { + int Err = 0; // default return code + + // Get the field metadata entries + std::shared_ptr AllMeta = Field::getFieldMetadata(FieldName); + + // Loop through all metadata - the Metadata type is an alias for a std::map + // so we can use map iterators for the loop + for (auto IMeta = AllMeta->begin(); IMeta != AllMeta->end(); ++IMeta) { + + // Get name + std::string MetaName = IMeta->first; + // Get value after determining the data type + std::any MetaVal = IMeta->second; + if (MetaVal.type() == typeid(I8)) { + I8 MetaValI8 = std::any_cast(MetaVal); + Err = IO::writeMeta(MetaName, MetaValI8, FileID, FieldID); + + } else if (MetaVal.type() == typeid(I4)) { + I4 MetaValI4 = std::any_cast(MetaVal); + Err = IO::writeMeta(MetaName, MetaValI4, FileID, FieldID); + + } else if (MetaVal.type() == typeid(R8)) { + R8 MetaValR8 = std::any_cast(MetaVal); + // if reduced precision is desired, convert to single (R4) + if (ReducePrecision) { + R4 MetaValR4 = MetaValR8; + Err = IO::writeMeta(MetaName, MetaValR4, FileID, FieldID); + } else { + Err = IO::writeMeta(MetaName, MetaValR8, FileID, FieldID); + } + + } else if (MetaVal.type() == typeid(R4)) { + R4 MetaValR4 = std::any_cast(MetaVal); + Err = IO::writeMeta(MetaName, MetaValR4, FileID, FieldID); + + } else if (MetaVal.type() == typeid(bool)) { + bool MetaValBool = std::any_cast(MetaVal); + Err = IO::writeMeta(MetaName, MetaValBool, FileID, FieldID); + + } else if (MetaVal.type() == typeid(std::string)) { + std::string MetaValStr = std::any_cast(MetaVal); + Err = IO::writeMeta(MetaName, MetaValStr, FileID, FieldID); + + // If the metadata was assigned using a string literal, std::any + // stores it as a char pointer so we need to convert it differently + } else if (MetaVal.type() == typeid(const char *)) { + const char *MetaValChar = std::any_cast(MetaVal); + Err = IO::writeMeta(MetaName, MetaValChar, FileID, FieldID); + + } else { // unknown data type + Err = 2; + LOG_ERROR("Unknown data type for Metadata {} in Field {}", MetaName, + FieldName); + return Err; + } + if (Err != 0) { + LOG_ERROR("Error trying to write Metadata {} for Field {}", MetaName, + FieldName); + return Err; + } + + } // end loop over metadata + + return Err; + +} // End writeFieldMeta + +//------------------------------------------------------------------------------ +// Write a field's data array, performing any manipulations to reduce +// precision or move data between host and device +int IOStream::writeFieldData( + std::shared_ptr FieldPtr, // [in] field to write + int FileID, // [in] id assigned to open file + int FieldID, // [in] id assigned to the field + std::map &AllDimIDs // [in] dimension IDs +) { + + int Err = 0; + + // Retrieve some basic field information + std::string FieldName = FieldPtr->getName(); + bool OnHost = FieldPtr->isOnHost(); + FieldType MyType = FieldPtr->getType(); + int NDims = FieldPtr->getNumDims(); + if (NDims < 1) { + LOG_ERROR("Invalid number of dimensions for Field {}", FieldName); + Err = 2; + return Err; + } + + // Create the decomposition needed for parallel I/O + int MyDecompID; + int LocSize; + std::vector DimLengths(NDims); + Err = computeDecomp(FieldPtr, AllDimIDs, MyDecompID, LocSize, DimLengths); + if (Err != 0) { + LOG_ERROR("Error computing decomposition for Field {}", FieldName); + return Err; + } + + // Extract and write the array of data based on the type, dimension and + // memory location. The IO routines require a contiguous data pointer on + // the host. Kokkos array types do not guarantee contigous memory for + // multi-dimensional arrays. Here we create a contiguous space and perform + // any other transformations (host-device data transfer, reduce precision). + void *DataPtr; + void *FillValPtr; + + // Vector for contiguous storage - the appropriate vector will be + // selected and resized later. + std::vector DataI4(1); + std::vector DataI8(1); + std::vector DataR4(1); + std::vector DataR8(1); + I4 FillValI4; + I8 FillValI8; + R4 FillValR4; + R8 FillValR8; + + switch (MyType) { + + // I4 Fields + case FieldType::I4: + + DataI4.resize(LocSize); + DataPtr = DataI4.data(); + // get fill value + Err = FieldPtr->getMetadata("FillValue", FillValI4); + if (Err != 0) { + LOG_ERROR("Error retrieving FillValue for Field {}", FieldName); + Err = 4; + return Err; + } + FillValPtr = &FillValI4; + + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DI4 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + DataI4[I] = Data(I); + } + } else { + Array1DI4 DataTmp = FieldPtr->getDataArray(); + HostArray1DI4 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + DataI4[I] = Data(I); + } + } + break; + case 2: + if (OnHost) { + HostArray2DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataI4[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } else { + Array2DI4 DataTmp = FieldPtr->getDataArray(); + HostArray2DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataI4[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } + break; + case 3: + if (OnHost) { + HostArray3DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataI4[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } else { + Array3DI4 DataTmp = FieldPtr->getDataArray(); + HostArray3DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataI4[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } + break; + case 4: + if (OnHost) { + HostArray4DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataI4[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } else { + Array4DI4 DataTmp = FieldPtr->getDataArray(); + HostArray4DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataI4[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } + break; + case 5: + if (OnHost) { + HostArray5DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataI4[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } else { + Array5DI4 DataTmp = FieldPtr->getDataArray(); + HostArray5DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataI4[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } + break; + + } // end switch NDims + break; // end I4 type + + // I8 Fields + case FieldType::I8: + + DataI8.resize(LocSize); + DataPtr = DataI8.data(); + // Get fill value + Err = FieldPtr->getMetadata("FillValue", FillValI8); + if (Err != 0) { + LOG_ERROR("Error retrieving FillValue for Field {}", FieldName); + Err = 4; + return Err; + } + FillValPtr = &FillValI8; + + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DI8 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + DataI8[I] = Data(I); + } + } else { + Array1DI8 DataTmp = FieldPtr->getDataArray(); + HostArray1DI8 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + DataI8[I] = Data(I); + } + } + break; + case 2: + if (OnHost) { + HostArray2DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataI8[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } else { + Array2DI8 DataTmp = FieldPtr->getDataArray(); + HostArray2DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataI8[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } + break; + case 3: + if (OnHost) { + HostArray3DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataI8[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } else { + Array3DI8 DataTmp = FieldPtr->getDataArray(); + HostArray3DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataI8[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } + break; + case 4: + if (OnHost) { + HostArray4DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataI8[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } else { + Array4DI8 DataTmp = FieldPtr->getDataArray(); + HostArray4DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataI8[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } + break; + case 5: + if (OnHost) { + HostArray5DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataI8[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } else { + Array5DI8 DataTmp = FieldPtr->getDataArray(); + HostArray5DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataI8[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } + break; + } // end switch NDims + break; // end I8 type + + // R4 Fields + case FieldType::R4: + + DataR4.resize(LocSize); + DataPtr = DataR4.data(); + // Get fill value + Err = FieldPtr->getMetadata("FillValue", FillValR4); + if (Err != 0) { + LOG_ERROR("Error retrieving FillValue for Field {}", FieldName); + Err = 4; + return Err; + } + FillValPtr = &FillValR4; + + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DR4 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + DataR4[I] = Data(I); + } + } else { + Array1DR4 DataTmp = FieldPtr->getDataArray(); + HostArray1DR4 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + DataR4[I] = Data(I); + } + } + break; + case 2: + if (OnHost) { + HostArray2DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataR4[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } else { + Array2DR4 DataTmp = FieldPtr->getDataArray(); + HostArray2DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataR4[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } + break; + case 3: + if (OnHost) { + HostArray3DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataR4[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } else { + Array3DR4 DataTmp = FieldPtr->getDataArray(); + HostArray3DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataR4[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } + break; + case 4: + if (OnHost) { + HostArray4DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataR4[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } else { + Array4DR4 DataTmp = FieldPtr->getDataArray(); + HostArray4DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataR4[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } + break; + case 5: + if (OnHost) { + HostArray5DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataR4[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } else { + Array5DR4 DataTmp = FieldPtr->getDataArray(); + HostArray5DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataR4[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } + break; + } // end switch NDims + break; // end R4 type + + // R8 Fields + case FieldType::R8: + + // Get fill value + Err = FieldPtr->getMetadata("FillValue", FillValR8); + if (Err != 0) { + LOG_ERROR("Error retrieving FillValue for Field {}", FieldName); + Err = 4; + return Err; + } + DataR8.resize(LocSize); + if (ReducePrecision) { + FillValR4 = FillValR8; + FillValPtr = &FillValR4; + DataR4.resize(LocSize); + DataPtr = DataR4.data(); + } else { + FillValPtr = &FillValR8; + DataPtr = DataR8.data(); + } + + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DR8 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + DataR8[I] = Data(I); + } + } else { + Array1DR8 DataTmp = FieldPtr->getDataArray(); + HostArray1DR8 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + DataR8[I] = Data(I); + } + } + break; + case 2: + if (OnHost) { + HostArray2DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataR8[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } else { + Array2DR8 DataTmp = FieldPtr->getDataArray(); + HostArray2DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + DataR8[VecAdd] = Data(J, I); + ++VecAdd; + } + } + } + break; + case 3: + if (OnHost) { + HostArray3DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataR8[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } else { + Array3DR8 DataTmp = FieldPtr->getDataArray(); + HostArray3DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + DataR8[VecAdd] = Data(K, J, I); + ++VecAdd; + } + } + } + } + break; + case 4: + if (OnHost) { + HostArray4DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataR8[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } else { + Array4DR8 DataTmp = FieldPtr->getDataArray(); + HostArray4DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + DataR8[VecAdd] = Data(L, K, J, I); + ++VecAdd; + } + } + } + } + } + break; + case 5: + if (OnHost) { + HostArray5DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataR8[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } else { + Array5DR8 DataTmp = FieldPtr->getDataArray(); + HostArray5DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + DataR8[VecAdd] = Data(M, L, K, J, I); + ++VecAdd; + } + } + } + } + } + } + break; + } // end switch NDims + if (ReducePrecision) { + for (int I = 0; I < LocSize; ++I) { + DataR4[I] = DataR8[I]; + } + } + break; // end R8 type + + default: + LOG_ERROR("Cannot determine data type for field {}", FieldName); + Err = 3; + break; + + } // end switch data type + + // Write the data + Err = OMEGA::IO::writeArray(DataPtr, LocSize, FillValPtr, FileID, MyDecompID, + FieldID); + if (Err != 0) { + LOG_ERROR("Error writing data array for field {} in stream {}", FieldName, + Name); + return Err; + } + + // Clean up the decomp + Err = OMEGA::IO::destroyDecomp(MyDecompID); + if (Err != 0) { + LOG_ERROR("Error destroying decomp for field {} in stream {}", FieldName, + Name); + return Err; + } + + return Err; + +} // end writeFieldData + +//------------------------------------------------------------------------------ +// Read a field's data array, performing any manipulations to reduce +// precision or move data between host and device +int IOStream::readFieldData( + std::shared_ptr FieldPtr, // [in] field to read + int FileID, // [in] id assigned to open file + std::map &AllDimIDs, // [in] dimension IDs + int &FieldID // [out] id assigned to the field +) { + + int Err = 0; + + // Retrieve some basic field information + std::string FieldName = FieldPtr->getName(); + // For MPAS back compatibility, the old name has a first letter that is + // lower case + std::string OldFieldName = FieldName; + OldFieldName[0] = std::tolower(OldFieldName[0]); + bool OnHost = FieldPtr->isOnHost(); + FieldType MyType = FieldPtr->getType(); + int NDims = FieldPtr->getNumDims(); + if (NDims < 1) { + LOG_ERROR("Invalid number of dimensions for Field {}", FieldName); + Err = 1; + return Err; + } + + // Compute the parallel decomposition + int DecompID; + int LocSize; + std::vector DimLengths(NDims); + Err = computeDecomp(FieldPtr, AllDimIDs, DecompID, LocSize, DimLengths); + if (Err != 0) { + LOG_ERROR("Error computing decomposition for Field {}", FieldName); + return Err; + } + + // The IO routines require a pointer to a contiguous memory on the host + // so we first read into a vector. Only one of the vectors below will + // be used and resized appropriately. + void *DataPtr; + std::vector DataI4(1); + std::vector DataI8(1); + std::vector DataR4(1); + std::vector DataR8(1); + + switch (MyType) { + case FieldType::I4: + DataI4.resize(LocSize); + DataPtr = DataI4.data(); + case FieldType::I8: + DataI8.resize(LocSize); + DataPtr = DataI8.data(); + case FieldType::R4: + DataR4.resize(LocSize); + DataPtr = DataR4.data(); + case FieldType::R8: + DataR8.resize(LocSize); + DataPtr = DataR8.data(); + } + + // read data into vector + Err = IO::readArray(DataPtr, LocSize, FieldName, FileID, DecompID, FieldID); + if (Err != 0) { + // For back compatibility, try to read again with old field name + Err = IO::readArray(DataPtr, LocSize, OldFieldName, FileID, DecompID, + FieldID); + if (Err == 0) { + LOG_INFO("Ignore PIO error for field {} ", FieldName); + LOG_INFO("Found field under old name {} ", OldFieldName); + } else { + LOG_ERROR("Error reading data array for {} in stream {}", FieldName, + Name); + return Err; + } + } + + // Unpack vector into array based on type, dims and location + switch (MyType) { + + // I4 Fields + case FieldType::I4: + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DI4 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataI4[I]; + } + } else { + Array1DI4 DataTmp = FieldPtr->getDataArray(); + HostArray1DI4 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataI4[I]; + } + deepCopy(DataTmp, Data); + } + break; + case 2: + if (OnHost) { + HostArray2DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + } else { + Array2DI4 DataTmp = FieldPtr->getDataArray(); + HostArray2DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + deepCopy(DataTmp, Data); + } + break; + case 3: + if (OnHost) { + HostArray3DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + } + } else { + Array3DI4 DataTmp = FieldPtr->getDataArray(); + HostArray3DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 4: + if (OnHost) { + HostArray4DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + } + } + } else { + Array4DI4 DataTmp = FieldPtr->getDataArray(); + HostArray4DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 5: + if (OnHost) { + HostArray5DI4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + } + } + } + } else { + Array5DI4 DataTmp = FieldPtr->getDataArray(); + HostArray5DI4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataI4[VecAdd]; + ++VecAdd; + } + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + } // end switch NDims + break; // end I4 fields + + // I8 Fields + case FieldType::I8: + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DI8 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataI8[I]; + } + } else { + Array1DI8 DataTmp = FieldPtr->getDataArray(); + HostArray1DI8 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataI8[I]; + } + deepCopy(DataTmp, Data); + } + break; + case 2: + if (OnHost) { + HostArray2DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + } else { + Array2DI8 DataTmp = FieldPtr->getDataArray(); + HostArray2DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + deepCopy(DataTmp, Data); + } + break; + case 3: + if (OnHost) { + HostArray3DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + } + } else { + Array3DI8 DataTmp = FieldPtr->getDataArray(); + HostArray3DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 4: + if (OnHost) { + HostArray4DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + } + } + } else { + Array4DI8 DataTmp = FieldPtr->getDataArray(); + HostArray4DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 5: + if (OnHost) { + HostArray5DI8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + } + } + } + } else { + Array5DI8 DataTmp = FieldPtr->getDataArray(); + HostArray5DI8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataI8[VecAdd]; + ++VecAdd; + } + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + } // end switch NDims + break; // end I8 fields + + // R4 Fields + case FieldType::R4: + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DR4 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataR4[I]; + } + } else { + Array1DR4 DataTmp = FieldPtr->getDataArray(); + HostArray1DR4 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataR4[I]; + } + deepCopy(DataTmp, Data); + } + break; + case 2: + if (OnHost) { + HostArray2DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + } else { + Array2DR4 DataTmp = FieldPtr->getDataArray(); + HostArray2DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + deepCopy(DataTmp, Data); + } + break; + case 3: + if (OnHost) { + HostArray3DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + } + } else { + Array3DR4 DataTmp = FieldPtr->getDataArray(); + HostArray3DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 4: + if (OnHost) { + HostArray4DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + } + } + } else { + Array4DR4 DataTmp = FieldPtr->getDataArray(); + HostArray4DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 5: + if (OnHost) { + HostArray5DR4 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + } + } + } + } else { + Array5DR4 DataTmp = FieldPtr->getDataArray(); + HostArray5DR4 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataR4[VecAdd]; + ++VecAdd; + } + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + } // end switch NDims + break; // end R4 fields + + // R8 Fields + case FieldType::R8: + switch (NDims) { + case 1: + if (OnHost) { + HostArray1DR8 Data = FieldPtr->getDataArray(); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataR8[I]; + } + } else { + Array1DR8 DataTmp = FieldPtr->getDataArray(); + HostArray1DR8 Data = createHostMirrorCopy(DataTmp); + for (int I = 0; I < DimLengths[0]; ++I) { + Data(I) = DataR8[I]; + } + deepCopy(DataTmp, Data); + } + break; + case 2: + if (OnHost) { + HostArray2DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + } else { + Array2DR8 DataTmp = FieldPtr->getDataArray(); + HostArray2DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int J = 0; J < DimLengths[0]; ++J) { + for (int I = 0; I < DimLengths[1]; ++I) { + Data(J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + deepCopy(DataTmp, Data); + } + break; + case 3: + if (OnHost) { + HostArray3DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + } + } else { + Array3DR8 DataTmp = FieldPtr->getDataArray(); + HostArray3DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int K = 0; K < DimLengths[0]; ++K) { + for (int J = 0; J < DimLengths[1]; ++J) { + for (int I = 0; I < DimLengths[2]; ++I) { + Data(K, J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 4: + if (OnHost) { + HostArray4DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + } + } + } else { + Array4DR8 DataTmp = FieldPtr->getDataArray(); + HostArray4DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int L = 0; L < DimLengths[0]; ++L) { + for (int K = 0; K < DimLengths[1]; ++K) { + for (int J = 0; J < DimLengths[2]; ++J) { + for (int I = 0; I < DimLengths[3]; ++I) { + Data(L, K, J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + case 5: + if (OnHost) { + HostArray5DR8 Data = FieldPtr->getDataArray(); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + } + } + } + } else { + Array5DR8 DataTmp = FieldPtr->getDataArray(); + HostArray5DR8 Data = createHostMirrorCopy(DataTmp); + int VecAdd = 0; + for (int M = 0; M < DimLengths[0]; ++M) { + for (int L = 0; L < DimLengths[1]; ++L) { + for (int K = 0; K < DimLengths[2]; ++K) { + for (int J = 0; J < DimLengths[3]; ++J) { + for (int I = 0; I < DimLengths[4]; ++I) { + Data(M, L, K, J, I) = DataR8[VecAdd]; + ++VecAdd; + } + } + } + } + } + deepCopy(DataTmp, Data); + } + break; + } // end switch NDims + break; // end R8 fields + + default: + LOG_ERROR("Invalid data type while reading field {} for stream {}", + FieldName, Name); + Err = 3; + return Err; + break; + + } // end switch data type + + // Clean up the decomp + Err = OMEGA::IO::destroyDecomp(DecompID); + if (Err != 0) { + LOG_ERROR("Error destroying decomp after reading field {} in stream {}", + FieldName, Name); + return Err; + } + + return Err; + +} // End readFieldData + +//------------------------------------------------------------------------------ +// Reads a stream if it is time. Returns an error code. This is the internal +// read function used by the public read interface. +int IOStream::readStream( + const Clock &ModelClock, // [in] model clock for getting time + Metadata &ReqMetadata, // [inout] global metadata to extract from file + bool ForceRead // [in] optional: read even if not time +) { + int Err = 0; // default return code + + // First check that this is an input stream + if (Mode != IO::ModeRead) { + LOG_ERROR("IOStream read: cannot read stream defined as output stream"); + Err = 1; + return Err; + } + + // If it is not time to read, return + if (!ForceRead) { + if (!MyAlarm.isRinging() and !OnStartup) + return Err; + if (UseStartEnd) { // If time outside interval, return + if (!StartAlarm.isRinging()) + return Err; + if (EndAlarm.isRinging()) + return Err; + } + } + + // Get current simulation time and time string + TimeInstant SimTime = ModelClock.getCurrentTime(); + std::string SimTimeStr = SimTime.getString(5, 0, "_"); + + // Reset alarms and flags + if (OnStartup) + OnStartup = false; + if (MyAlarm.isRinging()) + MyAlarm.reset(SimTime); + + // Create filename + std::string InFileName; + // If using pointer files for this stream, read the filename from the pointer + if (UsePointer) { + std::ifstream PtrFile(PtrFilename); + PtrFile >> InFileName; + PtrFile.close(); + } else if (FilenameIsTemplate) { + // create file name from template + InFileName = buildFilename(Filename, ModelClock); + } else { + InFileName = Filename; + } + + // Open input file + int InFileID; + Err = OMEGA::IO::openFile(InFileID, InFileName, Mode, IO::FmtDefault, + ExistAction); + if (Err != 0) { + LOG_ERROR("Error opening file {} for input", InFileName); + return Err; + } + + // Read any requested global metadata + for (auto Iter = ReqMetadata.begin(); Iter != ReqMetadata.end(); ++Iter) { + std::string MetaName = Iter->first; + std::any MetaTmp = Iter->second; + + I4 MetaValI4; + I8 MetaValI8; + R4 MetaValR4; + R8 MetaValR8; + bool MetaValBool; + std::string MetaValStr; + + // Read metadata based on type + int ErrRead = 0; + if (MetaTmp.type() == typeid(I8)) { + ErrRead = IO::readMeta(MetaName, MetaValI8, InFileID, IO::GlobalID); + ReqMetadata[MetaName] = MetaValI8; + } else if (MetaTmp.type() == typeid(I4)) { + ErrRead = IO::readMeta(MetaName, MetaValI4, InFileID, IO::GlobalID); + ReqMetadata[MetaName] = MetaValI4; + } else if (MetaTmp.type() == typeid(R8)) { + ErrRead = IO::readMeta(MetaName, MetaValR8, InFileID, IO::GlobalID); + ReqMetadata[MetaName] = MetaValR8; + } else if (MetaTmp.type() == typeid(R4)) { + ErrRead = IO::readMeta(MetaName, MetaValR4, InFileID, IO::GlobalID); + ReqMetadata[MetaName] = MetaValR4; + } else if (MetaTmp.type() == typeid(std::string)) { + ErrRead = IO::readMeta(MetaName, MetaValStr, InFileID, IO::GlobalID); + ReqMetadata[MetaName] = MetaValStr; + // If ReqMetadata was initialized with a string literal, we detect + // the type but replace it with a std::string + } else if (MetaTmp.type() == typeid(const char *)) { + ErrRead = IO::readMeta(MetaName, MetaValStr, InFileID, IO::GlobalID); + ReqMetadata[MetaName] = MetaValStr; + } else { + LOG_ERROR("Metadata read failed: unknown data type for {} in file {}", + MetaName, InFileName); + ErrRead = 2; + } + + if (ErrRead != 0) { + LOG_ERROR("Error reading metadata {} from file {}", InFileName); + return ErrRead; + } + } // end loop over requested metadata + + // Get dimensions from file and check that file has same dimension lengths + std::map AllDimIDs; + Err = defineAllDims(InFileID, AllDimIDs); + if (Err != 0) { + LOG_ERROR("Error defining dimensions for file {} ", InFileName); + return Err; + } + + // For each field in the contents, define field and read field data + for (auto IFld = Contents.begin(); IFld != Contents.end(); ++IFld) { + + // Retrieve the field name and pointer + std::string FieldName = *IFld; + std::shared_ptr ThisField = Field::get(FieldName); + + // Extract the data pointer and read the data array + int FieldID; // not currently used but available if field metadata needed + Err = readFieldData(ThisField, InFileID, AllDimIDs, FieldID); + if (Err != 0) { + LOG_ERROR("Error reading field data for Field {} in Stream {}", + FieldName, Name); + return Err; + } + + } // End loop over field list + + // Close input file + Err = IO::closeFile(InFileID); + if (Err != 0) { + LOG_ERROR("Error closing input file {}", InFileName); + return Err; + } + + LOG_INFO("Successfully read stream {} from file {}", Name, InFileName); + + // End of routine - return + return Err; + +} // End read + +//------------------------------------------------------------------------------ +// Writes stream. This is the internal member write function used by the +// public write interfaces. +int IOStream::writeStream( + const Clock &ModelClock, // [in] Model clock needed for time stamps + bool ForceWrite, // [in] Optional: write even if not time + bool FinalCall // [in] Optional flag if called from finalize +) { + + int Err = 0; // default return code + + // First check that this is an output stream + if (Mode != IO::ModeWrite) { + LOG_ERROR("IOStream write: cannot write stream defined as input stream"); + Err = 1; + return Err; + } + + // If it is not time to write, return + if (!ForceWrite) { + bool StartupShutdown = OnStartup or (OnShutdown and FinalCall); + if (!MyAlarm.isRinging() and !StartupShutdown) + return Err; + if (UseStartEnd) { // If time outside interval, return + if (!StartAlarm.isRinging()) + return Err; + if (EndAlarm.isRinging()) + return Err; + } + } + + // Get current simulation time and time string + TimeInstant SimTime = ModelClock.getCurrentTime(); + std::string SimTimeStr = SimTime.getString(4, 0, "_"); + + // Reset alarms and flags + if (OnStartup) + OnStartup = false; + if (MyAlarm.isRinging()) + MyAlarm.reset(SimTime); + + // Create filename + std::string OutFileName; + if (FilenameIsTemplate) { + // create file name from template + OutFileName = buildFilename(Filename, ModelClock); + } else { + OutFileName = Filename; + } + + // Open output file + int OutFileID; + Err = OMEGA::IO::openFile(OutFileID, OutFileName, Mode, IO::FmtDefault, + ExistAction); + if (Err != 0) { + LOG_ERROR("IOStream::write: error opening file {} for output", + OutFileName); + return Err; + } + + // Write Metadata for global metadata (Code and Simulation) + // Always add current simulation time to Simulation metadata + Err = writeFieldMeta(CodeMeta, OutFileID, IO::GlobalID); + if (Err != 0) { + LOG_ERROR("Error writing Code Metadata to file {}", OutFileName); + return Err; + } + std::shared_ptr SimField = Field::get(SimMeta); + // Add the simulation time - if it was added previously, remove and + // re-add the current time + if (SimField->hasMetadata("SimulationTime")) + Err = SimField->removeMetadata("SimulationTime"); + Err = SimField->addMetadata("SimulationTime", SimTimeStr); + if (Err != 0) { + LOG_ERROR("Error adding current sim time to output {}", OutFileName); + return Err; + } + Err = writeFieldMeta(SimMeta, OutFileID, IO::GlobalID); + if (Err != 0) { + LOG_ERROR("Error writing Simulation Metadata to file {}", OutFileName); + return Err; + } + + // Assign dimension IDs for all defined dimensions + std::map AllDimIDs; + Err = defineAllDims(OutFileID, AllDimIDs); + if (Err != 0) { + LOG_ERROR("Error defined dimensions for file {}", OutFileName); + return Err; + } + + // Define each field and write field metadata + std::map FieldIDs; + I4 NDims; + std::vector DimNames; + std::vector FieldDims; + for (auto IFld = Contents.begin(); IFld != Contents.end(); ++IFld) { + + // Retrieve the field pointer + std::string FieldName = *IFld; + std::shared_ptr ThisField = Field::get(FieldName); + + // Retrieve the dimensions for this field and determine dim IDs + NDims = ThisField->getNumDims(); + if (NDims < 1) { + LOG_ERROR("Invalid number of dimensions for Field {}", FieldName); + Err = 2; + return Err; + } + DimNames.resize(NDims); + FieldDims.resize(NDims); + Err = ThisField->getDimNames(DimNames); + if (Err != 0) { + LOG_ERROR("Error retrieving dimension names for Field {}", FieldName); + return Err; + } + for (int IDim = 0; IDim < NDims; ++IDim) { + std::string DimName = DimNames[IDim]; + FieldDims[IDim] = AllDimIDs[DimName]; + } + + // Determine the data type and convert to IODataType + // Reduce floating point precision if requested + IO::IODataType MyIOType = getFieldIOType(ThisField); + + // Define the field and assign a FieldID + int FieldID; + Err = defineVar(OutFileID, FieldName, MyIOType, NDims, FieldDims.data(), + FieldID); + if (Err != 0) { + LOG_ERROR("Error defining field {} in stream {}", FieldName, Name); + return Err; + } + FieldIDs[FieldName] = FieldID; + + // Now we can write the field metadata + Err = writeFieldMeta(FieldName, OutFileID, FieldID); + if (Err != 0) { + LOG_ERROR("Error writing field metadata for field {} in stream {}", + FieldName, Name); + return Err; + } + } + + // End define mode + Err = IO::endDefinePhase(OutFileID); + if (Err != 0) { + LOG_ERROR("Error ending define phase for stream {}", Name); + return Err; + } + + // Now write data arrays for all fields in contents + for (auto IFld = Contents.begin(); IFld != Contents.end(); ++IFld) { + + // Retrieve the field pointer and FieldID + std::string FieldName = *IFld; + std::shared_ptr ThisField = Field::get(FieldName); + int FieldID = FieldIDs[FieldName]; + + // Extract and write the data array + Err = this->writeFieldData(ThisField, OutFileID, FieldID, AllDimIDs); + if (Err != 0) { + LOG_ERROR("Error writing field data for Field {} in Stream {}", + FieldName, Name); + return Err; + } + } + + // Close output file + Err = IO::closeFile(OutFileID); + if (Err != 0) { + LOG_ERROR("Error closing output file {}", OutFileName); + return Err; + } + + // If using pointer files for this stream, write the filename to the pointer + // after the file is successfully written + if (UsePointer) { + std::ofstream PtrFile(PtrFilename, std::ios::trunc); + PtrFile << OutFileName << std::endl; + PtrFile.close(); + } + + LOG_INFO("Successfully wrote stream {} to file {}", Name, OutFileName); + + // End of routine - return + return Err; + +} // end writeStream + +//------------------------------------------------------------------------------ +// Removes a single IOStream from the list of all streams. +void IOStream::erase(const std::string &StreamName // Name of IOStream to remove +) { + AllStreams.erase(StreamName); // use the map erase function to remove +} // End erase + +//------------------------------------------------------------------------------ +// Private utility functions for read/write +//------------------------------------------------------------------------------ +//------------------------------------------------------------------------------ +// Determines the IO Data type to use for a given field, taking into +// account the field's type and any reduced precision conversion +IO::IODataType IOStream::getFieldIOType( + std::shared_ptr FieldPtr // [in] pointer to Field +) { + + IO::IODataType ReturnType; + + // Get the field data type + FieldType MyType = FieldPtr->getType(); + + // Determine IO data type based on field type and any reduced precision + // conversion + switch (MyType) { + case FieldType::I4: + ReturnType = IO::IOTypeI4; + break; + case FieldType::I8: + ReturnType = IO::IOTypeI8; + break; + case FieldType::R4: + ReturnType = IO::IOTypeR4; + break; + case FieldType::R8: + if (ReducePrecision) { + ReturnType = IO::IOTypeR4; + } else { + ReturnType = IO::IOTypeR8; + } + break; + default: + std::string FieldName = FieldPtr->getName(); + LOG_ERROR("Cannot determine data type for field {}", FieldName); + break; + } + + return ReturnType; + +} // end getIOType + +//------------------------------------------------------------------------------ +// Builds a filename based on time information and a filename template +// where special tokens are translated as: +// $SimTime = simulation time in form YYYY-MM-DD_hh.mm.ss +// $WallTime = actual wallclock time in form YYYY-MM-DD_hh.mm.ss +// and individual components of simulation time can be used to customize +// the time stamp using: +// $Y = year part of simulation time stamp +// $M = month part of simulation time stamp +// $D = day part of simulation time stamp +// $h = hour part of simulation time stamp +// $m = minute part of simulation time stamp +// $s = seconds part of simulation time stamp +std::string IOStream::buildFilename( + const std::string &FilenameTemplate, // [in] template string for name + const Clock &ModelClock // [in] model clock for sim time +) { + + // Start with input template + std::string Outfile = FilenameTemplate; + + // Check if wallclock time is requested in the template - check for + // multiple variations + size_t Pos = Outfile.find("$WallTime"); + if (Pos == std::string::npos) + Pos = Outfile.find("$Walltime"); + if (Pos == std::string::npos) + Pos = Outfile.find("$walltime"); + if (Pos == std::string::npos) + Pos = Outfile.find("$WALLTIME"); + // If wallclock time is requested, replace with wallclock string + // in format YYYY-MM-DD_HH.MM.SS + if (Pos != std::string::npos) { // wallclock string was found + // Get wallclock time and convert to string + std::time_t Walltime = std::time(nullptr); + std::tm *WalltimeInfo = std::localtime(&Walltime); + char Walltimestamp[20]; + std::strftime(Walltimestamp, sizeof(Walltimestamp), "%Y-%m-%d_%H.%M.%S", + WalltimeInfo); + std::string WalltimeStr = Walltimestamp; // convert to std::string + // Replace template string with wallclock string + Outfile.replace(Pos, 9, WalltimeStr); + } + + // For remaining options, we will need the simulation time + TimeInstant SimTime = ModelClock.getCurrentTime(); + + // For many template options we also need the components of the current + // sim time so we extract them here + I8 IYear; + I8 IMonth; + I8 IDay; + I8 IHour; + I8 IMin; + I8 ISecW; + I8 ISecN; + I8 ISecD; + int Err = SimTime.get(IYear, IMonth, IDay, IHour, IMin, ISecW, ISecN, ISecD); + if (Err != 0) { + LOG_ERROR("Error getting components of simulation time for filename"); + } + + // Convert these to strings (will pad with zeros if needed later) + std::string SYear = std::to_string(IYear); + std::string SMonth = std::to_string(IMonth); + std::string SDay = std::to_string(IDay); + std::string SHour = std::to_string(IHour); + std::string SMin = std::to_string(IMin); + std::string SSec = std::to_string(ISecW); + + // Now check for the $SimTime token and replace with a sim time string + Pos = Outfile.find("$SimTime"); + if (Pos == std::string::npos) + Pos = Outfile.find("$Simtime"); + if (Pos == std::string::npos) + Pos = Outfile.find("$simtime"); + if (Pos == std::string::npos) + Pos = Outfile.find("$SIMTIME"); + if (Pos != std::string::npos) { // Found the sim time token + // Convert the SimTime to a string + int YearLength = SYear.length(); + int MinWidth = 4; + int YearWidth = std::max(YearLength, MinWidth); + std::string SimTimeStr = SimTime.getString(YearWidth, 0, "_"); + Outfile.replace(Pos, 8, SimTimeStr); + } + + // Check for each of the other standard tokens and replace with + // appropriate strings with padding as necessary + int NPad = 0; + int SLength = 0; + + // Year is requested. The length is assumed to be at least 4 and + // padded as needed + if ((Pos = Outfile.find("$Y")) != std::string::npos) { + int SLength = SYear.length(); + if (SLength < 4) { + NPad = 4 - SLength; + SYear.insert(0, NPad, '0'); + } + Outfile.replace(Pos, 2, SYear); + } + + // Month is requested, pad to 2 digit month and insert into filename + if ((Pos = Outfile.find("$M")) != std::string::npos) { + SLength = SMonth.length(); + if (SLength < 2) { + NPad = 2 - SLength; + SMonth.insert(0, NPad, '0'); + } + Outfile.replace(Pos, 2, SMonth); + } + + // Day is requested, pad to 2 digit month and insert into filename + if ((Pos = Outfile.find("$D")) != std::string::npos) { + SLength = SDay.length(); + if (SLength < 2) { + NPad = 2 - SLength; + SDay.insert(0, NPad, '0'); + } + Outfile.replace(Pos, 2, SDay); + } + + // Hour is requested, pad to 2 digit hour and insert into filename + if ((Pos = Outfile.find("$h")) != std::string::npos) { + SLength = SHour.length(); + if (SLength < 2) { + NPad = 2 - SLength; + SHour.insert(0, NPad, '0'); + } + Outfile.replace(Pos, 2, SHour); + } + + // Minute is requested, pad to 2 digit minute and insert into filename + if ((Pos = Outfile.find("$m")) != std::string::npos) { + SLength = SMin.length(); + if (SLength < 2) { + NPad = 2 - SLength; + SMin.insert(0, NPad, '0'); + } + Outfile.replace(Pos, 2, SMin); + } + + // Second is requested. We have already rounded down to whole + // seconds, so pad that to two digits if needed + if ((Pos = Outfile.find("$s")) != std::string::npos) { + SLength = SSec.length(); + if (SLength < 2) { + NPad = 2 - SLength; + SSec.insert(0, NPad, '0'); + } + Outfile.replace(Pos, 2, SSec); + } + + // Add netcdf suffix + Outfile = Outfile + ".nc"; + + return Outfile; + +} // End buildFilename + +//------------------------------------------------------------------------------ +// Sets ReducePrecision flag based on an input string, performing string +// manipulation for case insensitive comparison + +void IOStream::setPrecisionFlag( + const std::string &PrecisionString ///< [in] precision from input YAML +) { + + // Set default value + ReducePrecision = false; + + // Convert input string to lowercase for easier comparison + std::string PrecComp = PrecisionString; + std::transform(PrecComp.begin(), PrecComp.end(), PrecComp.begin(), + [](unsigned char c) { return std::tolower(c); }); + + // Reduced precision if input is single + if (PrecComp == "single") { + ReducePrecision = true; + } else if (PrecComp == "double") { + ReducePrecision = false; + } else { + LOG_ERROR("Unknown precision {} in stream {}", PrecisionString, Name); + } + +} // End setPrecisionFlag + +} // namespace OMEGA +//===----------------------------------------------------------------------===// diff --git a/components/omega/src/infra/IOStream.h b/components/omega/src/infra/IOStream.h new file mode 100644 index 000000000000..5833ebdeb695 --- /dev/null +++ b/components/omega/src/infra/IOStream.h @@ -0,0 +1,257 @@ +#ifndef OMEGA_IOSTREAM_H +#define OMEGA_IOSTREAM_H +//===-- infra/IOStream.h - IO stream class ----------------------*- C++ -*-===// +// +/// \file +/// \brief Defines IO Stream class and methods +/// +/// This header defines classes and methods for IO Streams. IOStreams define +/// reading and writing of fields from/to a data file. Each stream contains +/// information on the file location, the frequency of input/output, the +/// fields to be read/written and other information necessary for file I/O. +/// Note that this is different from the C++ stdlib iostreams +/// +//===----------------------------------------------------------------------===// + +#include "Config.h" +#include "DataTypes.h" +#include "Dimension.h" +#include "Field.h" +#include "IO.h" +#include "Logging.h" +#include "TimeMgr.h" // need Alarms, TimeInstant +#include +#include +#include +#include + +namespace OMEGA { + +class IOStream { + + private: + /// Store and maintain all defined streams + static std::map> AllStreams; + + /// Private variables specific to a stream + std::string Name; ///< name of stream + std::string Filename; ///< filename or filename template (with path) + bool FilenameIsTemplate; ///< true if the filename is a template + IO::IfExists ExistAction; ///< action if file exists (write only) + + IO::Mode Mode; ///< mode (read or write) + bool ReducePrecision; ///< flag to use 32-bit precision for 64-bit floats + Alarm MyAlarm; ///< time mgr alarm for read/write + bool OnStartup; ///< flag to read/write on model startup + bool OnShutdown; ///< flag to read/write on model shutdown + + /// A pointer file is used if we wish OMEGA to read the name of the file + /// from another file. This is useful for writing the name of a restart + /// file to be picked up by the next job submitted so that the input + /// configuration does not need to change with each submission + bool UsePointer; ///< flag for using a pointer file + std::string PtrFilename; ///< name of pointer file + + /// Use a start and end time to define an interval in which stream is active + /// The start is inclusive but the end time is not. + bool UseStartEnd; ///< flag for using start, end times + Alarm StartAlarm; ///< alarm ringing if equal to or past start time + Alarm EndAlarm; ///< alarm ringing if equal to or past end time + + /// Contents of stream in the form of a set of Field names + std::set Contents; + + /// Flag to determine whether the Contents have been validated or not + bool Validated; + + //---- Private utility functions to support public interfaces + /// Creates a new stream and adds to the list of all streams, based on + /// options in the input model configuration. This routine is called by + /// the IOStreams initialize function. It requires an initialized model + /// clock so that stream alarm can be attached to this clock during creation. + static int create(const std::string &StreamName, ///< [in] name of stream + Config &StreamConfig, ///< [in] input stream configuration + Clock &ModelClock ///< [inout] Omega model clock + ); + + /// Define all dimensions used. Returns an error code as well as a map + /// of dimension names to defined dimension IDs. + int defineAllDims( + int FileID, ///< [in] id assigned to the IO file + std::map &AllDimIDs ///< [out] dim name, assigned ID + ); + + /// Computes the parallel decomposition (offsets) for a field. + /// Needed for parallel I/O + int computeDecomp( + std::shared_ptr FieldPtr, ///< [in] field + std::map &AllDimIDs, ///< [in] dimension IDs + int &DecompID, ///< [out] ID assigned to the defined decomposition + I4 &LocalSize, ///< [out] size of the local array for this field + std::vector &DimLengths // [out] local dim lengths + ); + + /// Private function that performs most of the stream read - called by the + /// public read method + int readStream( + const Clock &ModelClock, ///< [in] Model clock for alarms, time stamp + Metadata &ReqMetadata, ///< [inout] global metadata to extract from file + bool ForceRead = false ///< [in] Optional: read even if not time + ); + + /// Private function that performs most of the stream write - called by the + /// public write method + int writeStream( + const Clock &ModelClock, ///< [in] Model clock for alarms, time stamp + bool ForceWrite = false, ///< [in] Optional: write even if not time + bool FinalCall = false ///< [in] Optional flag for shutdown + ); + + /// Write all metadata associated with a field + int writeFieldMeta(std::string FieldName, ///< [in] metadata from field; + int FileID, ///< [in] id assigned to open file + int FieldID ///< [in] id assigned to the field + ); + + /// Write a field's data array, performing any manipulations to reduce + /// precision or move data between host and device + int + writeFieldData(std::shared_ptr FieldPtr, ///< [in] field to write + int FileID, ///< [in] id assigned to open file + int FieldID, ///< [in] id assigned to the field + std::map &AllDimIDs ///< [in] dimension IDs + ); + + /// Read a field's data array, performing any manipulations to reduce + /// precision or move data between host and device + int + readFieldData(std::shared_ptr FieldPtr, ///< [in] field to read + int FileID, ///< [in] id assigned to open file + std::map &AllDimIDs, ///< [in] dimension IDs + int &FieldID ///< [out] id assigned to the field + ); + + /// Determines the IO Data type to use for a given field, taking into + /// account the field's type and any reduced precision conversion + IO::IODataType getFieldIOType( + std::shared_ptr FieldPtr ///< [in] field to extract type + ); + + /// Builds a filename based on time information and a filename template + /// where special tokens are translated as: + /// $SimTime = simulation time in form YYYY-MM-DD_hh.mm.ss + /// $WallTime = actual wallclock time in form YYYY-MM-DD_hh.mm.ss + /// and individual components of simulation time can be used to customize + /// the time stamp using: + /// $Y = year part of simulation time stamp + /// $M = month part of simulation time stamp + /// $D = day part of simulation time stamp + /// $h = hour part of simulation time stamp + /// $m = minute part of simulation time stamp + /// $s = seconds part of simulation time stamp + static std::string buildFilename( + const std::string &FilenameTemplate, ///< [in] template string for name + const Clock &ModelClock ///< [in] model clock for sim time + ); + + /// Sets ReducePrecision flag based on an input string, performing string + /// manipulation for case insensitive comparison + void setPrecisionFlag( + const std::string &PrecisionString ///< [in] precision from input YAML + ); + + public: + //--------------------------------------------------------------------------- + /// Default empty constructor + IOStream(); + + //--------------------------------------------------------------------------- + /// Creates all streams defined in the input configuration file. This does + /// not validate the contents of the streams since the relevant Fields + /// may not have been defined yet. Returns an error code. + static int init(Clock &ModelClock ///< [inout] Omega model clock + ); + + //--------------------------------------------------------------------------- + /// Performs a final write of any streams that have the OnShutdown option and + /// then removes all streams to clean up. Returns an error code. + static int + finalize(const Clock &ModelClock ///< [in] Model clock needed for time stamps + ); + + //--------------------------------------------------------------------------- + /// Retrieves a previously defined stream by name. + static std::shared_ptr + get(const std::string &StreamName ///< [in] name of stream to retrieve + ); + + //--------------------------------------------------------------------------- + /// Adds a field to the contents of a stream. Because streams may be created + /// before all Fields have been defined, we only store the name. Validity + /// is either checked during read/write or can be checked using the validate + /// function. + void addField(const std::string &FieldName ///< [in] Name of field + ); + + //--------------------------------------------------------------------------- + /// Removes a field from the contents. Provided for symmetry, but not + /// often used. + void removeField(const std::string &FieldName ///< [in] Name of field + ); + + //--------------------------------------------------------------------------- + /// Validates an IOStream and its contents. If used, this must be called at + /// the end of initialization to ensure all Fields have been defined. + /// This function also replaces any group names by the list of field members + /// so that they can individually be checked and later read/write functions + /// have the complete list of fields. Returns true if all contents and + /// variables are valid. + bool validate(); + + //--------------------------------------------------------------------------- + /// Determines whether the contents of the stream have been validated + /// or not. + bool isValidated(); + + //--------------------------------------------------------------------------- + /// Validates all streams and their contents. If used, this must be called at + /// the end of initialization to ensure all Fields have been defined. + /// Returns true if all streams are valid. + static bool validateAll(); + + //--------------------------------------------------------------------------- + /// Reads a stream if it is time. Returns an error code. + static int read(const std::string &StreamName, ///< [in] Name of stream + const Clock &ModelClock, ///< [in] Model clock for time info + Metadata &ReqMetadata, ///< [inout] Metadata desired in file + bool ForceRead = false ///< [in] opt: read even if not time + ); + + //--------------------------------------------------------------------------- + /// Writes a stream if it is time. Returns an error code. + static int + write(const std::string &StreamName, ///< [in] Name of stream + const Clock &ModelClock, ///< [in] Model clock for time stamps + bool ForceWrite = false ///< [in] opt: write even if not time + ); + + //--------------------------------------------------------------------------- + /// Loops through all streams and writes them if it is time. This is + /// useful if most I/O is consolidated at one point (eg end of step). + static int + writeAll(const Clock &ModelClock ///< [in] Model clock for time stamps + ); + + //--------------------------------------------------------------------------- + /// Removes a single IOStream from the list of all streams. + /// That process also decrements the reference counters for the + /// shared pointers and removes them if those counters reach 0. + static void erase(const std::string &StreamName /// Name of stream to remove + ); + +}; // end class IOStream + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// +#endif // OMEGA_IOSTREAM_H diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index ff9b5bbedbdd..a11b3a65c3e6 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -234,6 +234,17 @@ add_omega_test( "-n;8" ) +################## +# IOStream test +################## + +add_omega_test( + IOSTREAM_TEST + testIOStream.exe + infra/IOStreamTest.cpp + "-n;8" +) + ##################### # TendencyTerms test ##################### diff --git a/components/omega/test/base/IOTest.cpp b/components/omega/test/base/IOTest.cpp index 10a10480f16f..ec2bf41674ec 100644 --- a/components/omega/test/base/IOTest.cpp +++ b/components/omega/test/base/IOTest.cpp @@ -345,6 +345,12 @@ int main(int argc, char *argv[]) { RetVal += 1; LOG_ERROR("IOTest: error writing global char metadata FAIL"); } + Err = OMEGA::IO::writeMeta("StringLiteral", "MyString", OutFileID, + OMEGA::IO::GlobalID); + if (Err != 0) { + RetVal += 1; + LOG_ERROR("IOTest: error writing char literal metadata FAIL"); + } // Define variables/arrays int VarIDCellI4; @@ -709,6 +715,19 @@ int main(int argc, char *argv[]) { LOG_INFO("IOTest: read/write file metadata string test FAIL"); } + std::string MyStringNew; + Err = OMEGA::IO::readMeta("StringLiteral", MyStringNew, InFileID, + OMEGA::IO::GlobalID); + if (Err != 0) { + RetVal += 1; + LOG_ERROR("IOTest: error reading file string literal FAIL"); + } + if (MyStringNew == "MyString") { + LOG_INFO("IOTest: read/write file metadata string literal test PASS"); + } else { + RetVal += 1; + LOG_INFO("IOTest: read/write file metadata string literal test FAIL"); + } // Read arrays OMEGA::HostArray2DI4 NewI4Cell("NewI4Cell", NCellsSize, NVertLevels); OMEGA::HostArray2DI8 NewI8Cell("NewI8Cell", NCellsSize, NVertLevels); diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp new file mode 100644 index 000000000000..3864c3ecefec --- /dev/null +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -0,0 +1,294 @@ +//===-- Test driver for OMEGA IO Streams -------------------------*- C++ -*-===/ +// +/// \file +/// \brief Test driver for OMEGA IO Streams +/// +/// This driver tests the ability to write to/from files within Omega. An +/// IOStream is defined for each unique combination of read/write frequency, +/// contents of the file and other properties. +// +//===-----------------------------------------------------------------------===/ + +#include "IOStream.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "Field.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "TimeMgr.h" +#include "TimeStepper.h" +#include "mpi.h" +#include + +using namespace OMEGA; + +//------------------------------------------------------------------------------ +// Set some constant reference values for simplicity +const I4 RefI4 = 3; +const I8 RefI8 = 400000000; +const R4 RefR4 = 5.1; +const R8 RefR8 = 6.123456789; +const std::string RefStr = "Reference String"; + +//------------------------------------------------------------------------------ +// A simple test evaluation function +template +void TestEval(const std::string &TestName, T TestVal, T ExpectVal, int &Error) { + + if (TestVal == ExpectVal) { + LOG_INFO("{}: PASS", TestName); + } else { + LOG_ERROR("{}: FAIL", TestName); + ++Error; + } +} +//------------------------------------------------------------------------------ +// Initialization routine to create reference Fields +int initIOStreamTest(std::shared_ptr &ModelClock, // Model clock + Calendar &ModelCalendar) { + + int Err = 0; + int Err1 = 0; + int ErrRef = 0; + + // Initialize maching environment and logging + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + initLogging(DefEnv); + + // Read the model configuration + Config Config("omega"); + Err1 = Config::readAll("omega.yml"); + TestEval("Config read all", Err1, ErrRef, Err); + + // Initialize base-level IO + Err1 = IO::init(DefComm); + TestEval("IO Initialization", Err1, ErrRef, Err); + + // Initialize decomposition + Decomp::init(); + Decomp *DefDecomp = Decomp::getDefault(); + + // Initialize Halo updates + Halo::init(); + OMEGA::Halo *DefHalo = OMEGA::Halo::getDefault(); + + // Initialize Field + Err1 = Field::init(); + TestEval("IO Field initialization", Err1, ErrRef, Err); + + // Create the model clock and time step + TimeInstant SimStartTime(&ModelCalendar, 0001, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(2, TimeUnits::Hours); + ModelClock = std::make_shared(SimStartTime, TimeStep); + + // Initialize IOStreams + Err1 = IOStream::init(*ModelClock); + TestEval("IOStream Initialization", Err1, ErrRef, Err); + + // Initialize HorzMesh - this should read Mesh stream + Err1 = HorzMesh::init(); + TestEval("Horizontal mesh initialization", Err1, ErrRef, Err); + HorzMesh *DefMesh = HorzMesh::getDefault(); + I4 NCellsSize = DefMesh->NCellsSize; + + // Set vertical levels and time levels + I4 NVertLevels = 60; + std::shared_ptr VertDim = + Dimension::create("NVertLevels", NVertLevels); + + // Initialize time stepper needed before ocean state (for time levels) + Err1 = TimeStepper::init(); + TestEval("Ocean time step initialization", Err1, ErrRef, Err); + + // Initialize State + Err1 = OceanState::init(); + TestEval("Ocean state initialization", Err1, ErrRef, Err); + + // Add some global (Model and Simulation) metadata + std::shared_ptr CodeField = Field::get(CodeMeta); + std::shared_ptr SimField = Field::get(SimMeta); + + Err1 = CodeField->addMetadata("CodeIntTest", 3); + TestEval("Add code metadata int", Err1, ErrRef, Err); + Err1 = CodeField->addMetadata("CodeRealTest", 4.567); + TestEval("Add code metadata real", Err1, ErrRef, Err); + Err1 = CodeField->addMetadata("CodeBoolTest", true); + TestEval("Add code metadata bool", Err1, ErrRef, Err); + std::string CodeStrVal = "ASampleString"; + Err1 = CodeField->addMetadata("CodeStrTest", CodeStrVal); + TestEval("Add code metadata str", Err1, ErrRef, Err); + Err1 = CodeField->addMetadata("CodeVersion", "V0.0"); + TestEval("Add code metadata str literal", Err1, ErrRef, Err); + Err1 = SimField->addMetadata("ExpName", "IOStreamsTest"); + TestEval("Add ExpName metadata", Err1, ErrRef, Err); + std::string StartTimeStr = SimStartTime.getString(4, 2, "_"); + Err1 = SimField->addMetadata("SimStartTime", StartTimeStr); + TestEval("Add SimStartTime metadata", Err1, ErrRef, Err); + + // Define temperature and salinity tracer fields and create a tracer + // group + + std::vector DimNames(2); + DimNames[0] = "NCells"; + DimNames[1] = "NVertLevels"; + + // 2D Fields on device + + DimNames[0] = "NCells"; + DimNames[1] = "NVertLevels"; + Real FillValue = -1.2345e-30; + auto TempField = Field::create( + "Temperature", "Potential temperature at cell centers", "deg C", + "sea_water_pot_tem", -3.0, 100.0, FillValue, 2, DimNames); + auto SaltField = + Field::create("Salinity", "Salinity at cell centers", "", + "sea_water_salinity", 0.0, 100.0, FillValue, 2, DimNames); + + // Create Tracer group + auto TracerGroup = FieldGroup::create("Tracers"); + + // Add fields to tracer group + Err1 = TracerGroup->addField("Temperature"); + TestEval("Add Temperature to tracer group", Err1, ErrRef, Err); + Err1 = TracerGroup->addField("Salinity"); + TestEval("Add Salinity to tracer group", Err1, ErrRef, Err); + + // Also create a Restart group + auto RestartGroup = FieldGroup::create("Restart"); + + // Add fields to restart group + Err1 = RestartGroup->addField("Temperature"); + TestEval("Add Temperature to restart group", Err1, ErrRef, Err); + Err1 = RestartGroup->addField("Salinity"); + TestEval("Add Salinity to restart group", Err1, ErrRef, Err); + + // End of init + return Err; + +} // End initialization IOStream test + +//------------------------------------------------------------------------------ +// We will test the IOStream interfaces by defining Fields and reading stream +// configurations during init. We will test periodic writing of data and +// reading/writing restart and initial data. + +int main(int argc, char **argv) { + + int Err = 0; + int Err1 = 0; + int ErrRef = 0; + + // Initialize the global MPI and Kokkos environments + MPI_Init(&argc, &argv); + Kokkos::initialize(); + { + + Calendar CalGreg("Gregorian", OMEGA::CalendarGregorian); + std::shared_ptr ModelClock; + + // Call initialization to create reference IO field + Err1 = initIOStreamTest(ModelClock, CalGreg); + TestEval("Initialize IOStream test", Err1, ErrRef, Err); + + // Retrieve dimension lengths and some mesh info + I4 NCellsSize = Dimension::getDimLengthLocal("NCells"); + I4 NVertLevels = Dimension::getDimLengthLocal("NVertLevels"); + Decomp *DefDecomp = Decomp::getDefault(); + I4 NCellsOwned = DefDecomp->NCellsOwned; + Array1DI4 CellID = DefDecomp->CellID; + + // Create data arrays + + Array2DR8 Temp("Temp", NCellsSize, NVertLevels); + Array2DR8 Salt("Salt", NCellsSize, NVertLevels); + Array2DR8 Test("Test", NCellsSize, NVertLevels); + + // Attach data arrays to fields + + Err1 = Field::attachFieldData("Temperature", Temp); + TestEval("Attach temperature data to field", Err1, ErrRef, Err); + Err1 = Field::attachFieldData("Salinity", Salt); + TestEval("Attach salinity data to field", Err1, ErrRef, Err); + + // Validate all streams (Mesh stream already validated in HorzMesh?) + bool AllValidated = IOStream::validateAll(); + TestEval("IOStream Validation", AllValidated, true, Err); + + // Read restart file for initial temperature and salinity data + Metadata ReqMetadata; // leave empty for now - no required metadata + Err1 = IOStream::read("InitialState", *ModelClock, ReqMetadata); + TestEval("Read restart file", Err1, ErrRef, Err); + + // Overwrite salinity array with values associated with global cell + // ID to test proper indexing of IO + parallelFor( + {NCellsSize, NVertLevels}, KOKKOS_LAMBDA(int Cell, int K) { + Salt(Cell, K) = 0.0001_Real * (CellID(Cell) + K); + Test(Cell, K) = Salt(Cell, K); + }); + + // Create a stop alarm at 1 year for time stepping + TimeInstant StopTime(&CalGreg, 0002, 1, 1, 0, 0, 0.0); + Alarm StopAlarm("Stop Time", StopTime); + Err1 = ModelClock->attachAlarm(&StopAlarm); + TestEval("Attach stop alarm", Err1, ErrRef, Err); + + // Overwrite + // Step forward in time and write files if it is time + while (!StopAlarm.isRinging()) { + ModelClock->advance(); + TimeInstant CurTime = ModelClock->getCurrentTime(); + std::string CurTimeStr = CurTime.getString(4, 2, " "); + + Err1 = IOStream::writeAll(*ModelClock); + if (Err1 != 0) // to prevent too much output in log + TestEval("Write all streams " + CurTimeStr, Err1, ErrRef, Err); + } + + // Force read the latest restart and check the results + // We use this log message to act as a barrier/delay to make sure the + // restart write has finished before we read. + LOG_INFO("Test reading final restart"); + bool ForceRead = true; + Err1 = IOStream::read("RestartRead", *ModelClock, ReqMetadata, ForceRead); + TestEval("Restart force read", Err1, ErrRef, Err); + + Err1 = 0; + auto DataReducer = Kokkos::Sum(Err1); + + parallelReduce( + {NCellsOwned, NVertLevels}, + KOKKOS_LAMBDA(int Cell, int K, I4 &Err1) { + if (Salt(Cell, K) != Test(Cell, K)) + ++Err1; + }, + DataReducer); + TestEval("Check Salt array ", Err1, ErrRef, Err); + + // Write final output and remove all streams + IOStream::finalize(*ModelClock); + } + + // Clean up environments + OceanState::clear(); + HorzMesh::clear(); + Field::clear(); + Dimension::clear(); + Decomp::clear(); + Kokkos::finalize(); + MPI_Finalize(); + + if (Err >= 256) + Err = 255; + + // End of testing + return Err; +} +//===--- End test driver for IO Streams ------------------------------------===/