From 2e30794c14b92c11c794fcad88598a121e80ead1 Mon Sep 17 00:00:00 2001 From: Youngsung Kim Date: Wed, 18 Sep 2024 13:08:39 -0400 Subject: [PATCH] Add Tracers infrastructure Created the Tracers class and implemented its methods. Unit tests and documentation were also added. --- components/omega/configs/Default.yml | 3 + components/omega/doc/design/Tracers.md | 566 +++++++++++----------- components/omega/doc/devGuide/Tracers.md | 210 ++++++++ components/omega/doc/index.md | 2 + components/omega/doc/userGuide/Tracers.md | 42 ++ components/omega/src/ocn/TracerDefs.inc | 25 + components/omega/src/ocn/Tracers.cpp | 527 ++++++++++++++++++++ components/omega/src/ocn/Tracers.h | 212 ++++++++ components/omega/test/CMakeLists.txt | 11 + components/omega/test/ocn/TracersTest.cpp | 491 +++++++++++++++++++ 10 files changed, 1803 insertions(+), 286 deletions(-) create mode 100644 components/omega/doc/devGuide/Tracers.md create mode 100644 components/omega/doc/userGuide/Tracers.md create mode 100644 components/omega/src/ocn/TracerDefs.inc create mode 100644 components/omega/src/ocn/Tracers.cpp create mode 100644 components/omega/src/ocn/Tracers.h create mode 100644 components/omega/test/ocn/TracersTest.cpp diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 2fb2803081d1..8c03bc156857 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -25,3 +25,6 @@ Omega: ViscDel2: 1.0e3 VelHyperDiffTendencyEnable: true ViscDel4: 1.2e11 + Tracers: + Base: [Temp, Salt] + Debug: [Debug1, Debug2, Debug3] diff --git a/components/omega/doc/design/Tracers.md b/components/omega/doc/design/Tracers.md index 82ba3377ffa1..827d6771f3fa 100644 --- a/components/omega/doc/design/Tracers.md +++ b/components/omega/doc/design/Tracers.md @@ -1,415 +1,409 @@ (omega-design-tracers)= -# Tracer infrastructure + +# Tracer Infrastructure ## 1 Overview -Tracers are either heat (temperature) or material carried by a -fluid parcel (eg salt, chemical or biological constituents). Here -we describe a tracer module that manages how tracers are defined, -described, stored and accessed throughout OMEGA. For reasons -outlined in the requirements below, we do not address -the specific methods or algorithms for transporting or mixing -tracers here. Those algorithms will be covered in other +Tracers refer to either heat (temperature) or material carried by a fluid +parcel (e.g., salt, chemical, or biological constituents). This document +describes a tracer module that manages the definition, description, storage, +and access of tracers throughout OMEGA. For reasons outlined in the following +requirements, specific methods or algorithms for transporting or mixing +tracers are not addressed here. These algorithms will be covered in other design documents. ## 2 Requirements ### 2.1 Requirement: Tracer definition and metadata -The tracer module must provide a single location in which all -supported tracers are defined. This definition includes all -metadata associated with the tracer - see related documents -for metadata and IO. Where possible, CF compliant metadata -must be used. +The tracer module must provide a centralized location where all supported +tracers are defined. This definition includes all metadata associated with +the tracers—see related documents for metadata and I/O. Where possible, +CF-compliant metadata must be used. ### 2.2 Requirement: Tracer identification -The tracer module must provide a simple way to identify a tracer, -either by name or by index. The latter may be preferred for -performance reasons. +The tracer module must provide a simple way to identify a tracer, either +by name or by index. The latter may be preferred for performance reasons. ### 2.3 Requirement: Tracer groups -It is often convenient to select or define properties of tracers -as a group. It must be possible to assign tracers to a particular -group or groups. Properties of the group should apply to all -tracers assigned to that group. Tracers can be members of multiple -groups as long as properties of each group do not conflict. There -will be a tracer group for All that will include all selected -tracers. For performance reasons, it may be desirable for some -non-intersecting groups to be assigned contiguous memory or with -a contiguous index space. +It is often convenient to select or define properties for tracers as +a group. It must be possible to assign tracers to particular groups. +Properties of the group should apply to all tracers assigned to that group. +Tracers can belong to multiple groups, as long as properties of each group +do not conflict. A group called "All" will include all selected tracers. +For performance reasons, it may be desirable for some non-intersecting +groups to be assigned contiguous memory or contiguous index space. ### 2.4 Requirement: Tracer selection -At the start of the simulation, users must be able to select -which of the supported tracers are enabled for the experiment -through the configuration input file. For a more concise -configuration, it should be possible to select tracers by group -as well as individually. +At the start of a simulation, users must be able to select which of the +supported tracers are enabled for the experiment via the configuration +input file. To simplify the configuration, tracers should be selectable +by group as well as individually. -### 2.5 Requirement: Tracer restart and IO +### 2.5 Requirement: Tracer restart and I/O -With few exceptions, tracers will generally be part of all -restarts and will often be included in other output for -postprocessing and analysis. Practically, this is a -requirement for tracers to be registered as available IO -fields (see related IO designs) so that they can be read -from or written to a file. +With few exceptions, tracers will generally be included in all restarts and +often in other outputs for post-processing and analysis. Thus, it is +essential for tracers to be registered as available I/O fields (see related +I/O designs) so that they can be read from or written to a file. ### 2.6 Requirement: Time levels -Time integration schemes often require the storage of multiple -time levels. Tracer storage must support multiple time levels and -a mechanism for specifying which time level is needed for a -given calculation. Ideally the change of time levels at the end -of a time step should avoid a memory copy and be a simple change -of time index or change of pointer. +Time integration schemes often require the storage of multiple time levels. +Tracer storage must support multiple time levels, and a mechanism for +specifying which time level is needed for a given calculation must be +provided. Ideally, changing time levels at the end of a time step should +avoid a memory copy and instead involve a simple change of time index or +pointer. ### 2.7 Requirement: Acceleration or supercycling -Some tracers or tracer groups can be evolved using a longer -time step. For performance, it should be possible to advance -these tracers with a longer time step and avoid computing -tendencies at every dynamical time step. - -### 2.7 Desired: Per-tracer/group algorithmic requirements +Some tracers or tracer groups can be evolved using a longer time step. For +performance reasons, it should be possible to advance these tracers with a +longer time step and avoid computing tendencies at every dynamical time step. -Different tracers may require different methods for their evolution -to preserve desired properties. For example, temperature and -salinity should preserve tracer consistency with the continuity -equation. Many tracers will require monotone or positive-definite -advection schemes to prevent negative concentrations. We might -even use a different time stepping scheme for supercycled tracers. +### 2.8 Desired: Per-tracer/group algorithmic requirements +Different tracers may require different methods for their evolution to +preserve desired properties. For example, temperature and salinity should +maintain consistency with the continuity equation. Many tracers will require +monotonic or positive-definite advection schemes to prevent negative +concentrations. We may even use a different time-stepping scheme for +supercycled tracers. ## 3 Algorithmic Formulation -Algorithms used to advect, mix and evolve tracers in time are -described elsewhere in the respective modules that compute these -tendencies. +The algorithms used to advect, mix, and evolve tracers in time are described +elsewhere in the respective modules responsible for computing these tendencies. ## 4 Design -We will store all tracers in a single static array on the device: +We will store all tracers in a vector of 3-dimensional device arrays. + ```c++ - Array4DReal AllTracers(MaxTimeLvls, NumTracers, NCellsSize, - VertSize); + std::vector TracerArrays; // multiple time levels ``` -This will allow us to parallelize over tracer and cell dimensions -with the vertical dimension either parallel or vectorized depending -on architecture and computation involved. Many of the requirements -for grouping, selecting and swapping time levels can be achieved -through collection and manipulation of array indices. Individual -tracers can still be extracted as arrays with contiguous memory. - -All tracers will be defined in a set of code blocks (described -below) in a file `TracerDefs.inc` that will be directly included -in the init routine to define all supported tracers and their -metadata. Tracer groups are defined as a set of tracer indices -that belong to the group. A set of primary groups will be defined -in `TracerDefs.inc` that define membership for the purposes of -selecting tracers. Additional groups may be added by other -modules as needed (eg if we decide to do per-tracer algorithm -choices, the tracer advection module may define groups that -assign tracers to a particular method). - -### 4.1 Data types and parameters + +Each device array in the vector has dimensions corresponding to the number +of tracers, the number of all cells in the rank, and the vertical length. +This allows parallelization over tracer and cell dimensions, with the vertical +dimension either parallel or vectorized, depending on the architecture and +computation involved. Many of the requirements for grouping, selecting, and +swapping time levels can be achieved through the collection and manipulation +of array indices. Individual tracers can still be extracted as arrays with +contiguous memory. + +All tracers will be defined in a set of code blocks (described below) in +a file named `TracerDefs.inc`, which will be directly included in the +initialization routine to define all supported tracers and their metadata. +Tracer groups are defined as sets of tracer indices in the OMEGA YAML +configuration file. + +### 4.1 Data Types and Parameters #### 4.1.1 Definition -All supported tracers are defined in a file called `TracerDefs.inc` -with the code blocks needed to define the tracer and its metadata. -This must include all tracers potentially used in the model. -The file will look something like: +All supported tracers are defined in a file called `TracerDefs.inc`, which +includes the code blocks needed to define each tracer and its metadata. This +file must include all tracers potentially used in the model. The file will +look something like this: + ```c++ -// At the top of the file, list all available tracers in a -// comment to provide a sort of table of contents. +// At the top of the file, list all available tracers in a comment +// to provide a sort of table of contents. // Tracers: // Temp // Salt // [Others...] -// Tracer Groups: -// Base (or Active?) -// Debug -// Ecosys -// Analysis - -// Individual tracer entries will look like: -// Tracer: Temp - - Err = OMEGA::Tracer::define( - "Temp", // Name of var - "Potential Temperature", // long Name or description - "degree_C", // units - "sea_water_potential_temperature", /// CF standard Name - -273.15, // min valid field value - 100.0, // max valid field value - 1.e33 // fill value for undef entries - ); - Err = OMEGA::Tracer::addToGroup("Base", "Temp"); + +define( + "Temp", // Name of variable + "Potential Temperature", // Long name or description + "degree_C", // Units + "sea_water_potential_temperature", // CF standard name + -273.15, // Min valid field value + 100.0, // Max valid field value + 1.e33 // Fill value for undefined entries +); ``` #### 4.1.2 Configuration -Users will select desired tracers either individually or by -group. If the group is in the list, it implies the entire -group, but the user can choose a subset of a group by -supplying the individual tracer names and omitting the group name. +Users will select desired tracers either individually or by group. If a group +is listed, it implies the entire group, but users can choose a subset of +a group by specifying individual tracer names and omitting the group name. + ```yaml omega: Tracers: - - Base (or Active or ? - basically the required T,S) - - Debug - - Ecosys - - IdealAge - - [other individuals or groups as needed] - + Base: [Temp, Salt] + Debug: [Debug1, Debug2, Debug3] + [other individual tracers or groups as needed] ``` -For requirement 2.7 (per-tracer algorithm choice), the configuration -could look something like: +For requirement 2.7 (per-tracer algorithm choice), the configuration might +look something like this: + ```yaml omega: TracerAdv: - # start with default values for all tracers - - All + # Start with default values for all tracers + - All: Method: standard SomeParameter: 0.5 - # specialized options - override default by group or tracer name + # Specialized options - override default by group or tracer name - Ecosys: Method: MyMonotoneMethod - - [Other tracers or tracer groups, only if needed] - Method: My additional option + - [Other tracers or tracer groups, only if needed]: + Method: MyAdditionalOption SomeParameter: 1.0 ``` +#### 4.1.3 Classes/Structs/Data Types -#### 4.1.3 Class/structs/data types - -There is a single class that contains the static array for storing -all tracers and a number of containers to manage tracer groups and -indices. Methods are described in section 4.2 below. +There will be a single class that contains the static array for storing +all tracers, along with several containers for managing tracer groups and +indices. Methods are described in Section 4.2 below. ```c++ +/// The Tracers class provides a container for tracer arrays and methods class Tracers { - private: - // Total number of tracers activated for this simulation - static int NumTracers - - // static storage of the tracer array for all tracers - static Array4DReal AllTracers; - - // maps for matching tracer names with indices (both directions) - static std::map Index - static std::map Name - - // maps for managing tracer groups - /// This map associates a group name with a list of tracer - /// indices that are assigned to the group - static std::map> Groups - - /// This map matches the group name with a boolean vector - /// of NumTracers that flags a tracer index as a member - static std::map> MemberFlag; + private: + + // Total number of tracers defined at initialization + static I4 NumTracers; + + // Static storage of the tracer arrays for device and host memory spaces + // The arrays are 3-dimensional, with dimensions for Tracer, Cell, and + // Vertical layers + static std::vector + TracerArrays; ///< Time levels -> [Tracer, Cell, Vertical] + static std::vector + TracerArraysH; ///< Time levels -> [Tracer, Cell, Vertical] + + // Maps for managing tracer groups + // The key of this map is a group name and the value is a pair of + // GroupStartIndex and GroupLength + static std::map> TracerGroups; + + // Maps for matching tracer names with indices (in both directions) + static std::map TracerIndexes; + static std::map TracerNames; + + // Halo exchange + static Halo *MeshHalo; + + // Tracer dimension names + static std::vector TracerDimNames; + + // Current time index + // This index is circular, so it returns to index 0 when it exceeds the max + // index + static I4 CurTimeIndex; ///< Time dimension array index for the current level + + // Pack tracer field name + static std::string packTracerFieldName(const std::string &TracerName); + + // Defines all tracers locally without allocating memory + static I4 + define(const std::string &Name, ///< [in] Name of the tracer + const std::string &Description, ///< [in] Long name or description + const std::string &Units, ///< [in] Units + const std::string &StdName, ///< [in] CF standard name + const Real ValidMin, ///< [in] Min valid field value + const Real ValidMax, ///< [in] Max valid field value + const Real FillValue ///< [in] Value for undefined entries + ); public: [methods defined below] -} +}; ``` +Here is the corrected version with improved grammar and clarity: + +--- ### 4.2 Methods -#### 4.2.1 init +#### 4.2.1 `init` + +The initialization method defines all the tracers by including the +`TracerDefs.inc` file. Once all fields are defined, it allocates the full +tracer array based on the input configuration of selected tracers. The +interface is as follows: -The initialization method is where all the tracers are defined -by including the `TracerDefs.inc` file. Once all the fields have -been defined, it allocates the full tracer array based on the -input configuration of selected tracers. The interface is simply: ```c++ static int Err = init(); ``` -Note that this function does not actually initialize the -tracer values. The values will be filled either by input IO streams -or by individual modules associated with some tracers and tracer -groups. - -#### 4.2.2 define - -All tracers must be defined. This routine takes as input the -name of the tracer and some required metadata. It uses these -inputs to create a MetaData entry and an IOField entry, but only -if the tracer is selected in the input configuration file. It also -increments the number of tracers and assigns a tracer index if the -tracer is selected. It returns an error code that is non-zero if -any of a number of error conditions are encountered. + +Note that this function does not initialize the tracer values. The values +will be filled either by input I/O streams or by individual modules +associated with specific tracers and tracer groups. + +#### 4.2.2 `define` + +All tracers must be defined. This routine takes the name of the tracer and +some required metadata as input. It uses this information to create +a `MetaData` entry and an `IOField` entry, but only if the tracer is +selected in the input configuration file. It also increments the number of +tracers and assigns a tracer index if the tracer is selected. The function +returns a non-zero error code if any of several possible error conditions +are encountered. + ```c++ static int define( - const std::string &Name, ///< [in] Name of tracer + const std::string &Name, ///< [in] Name of the tracer const std::string &Description, ///< [in] Long name or description const std::string &Units, ///< [in] Units - const std::string &StdName, ///< [in] CF standard Name - OMEGA::Real ValidMin, ///< [in] min valid field value - OMEGA::Real ValidMax, ///< [in] max valid field value - OMEGA::Real FillValue ///< [in] value for undef entries + const std::string &StdName, ///< [in] CF standard name + OMEGA::Real ValidMin, ///< [in] Minimum valid field value + OMEGA::Real ValidMax, ///< [in] Maximum valid field value + OMEGA::Real FillValue ///< [in] Fill value for undefined entries ); ``` -#### 4.2.3 addToGroup +#### 4.2.3 `getNumTracers` + +This function retrieves the total number of tracers activated for the current +simulation. -Tracers can be added to a group with this method: -```c++ -static int addToGroup( - const std::string &GroupName, ///< [in] name of group - const std::string &TracerName ///< [in] name of tracer to add -); -``` -and this routine is typically called right after the tracer is -defined and included as part of the `TracerDefs.inc` file. Tracers -can be added to more than one group so multiple calls are possible. -If the group does not yet exist, one will be created with the name -provided. If the tracer already belongs to the group, no action is -taken. At this point, not all information on tracers and groups are -available since other fields in a group may not have been defined -or added yet. Once all fields have been defined and added to groups -(i.e. at the end of the tracers init routine), the `pruneGroups` -function is called to ensure the groups and group members are -consistent with the tracers selected in the input config file. - -#### 4.2.4 organize - -During the definition phase, groups are defined that may not -have been fully selected and some information cannot be inferred. -A final pass of the tracers and groups will be necessary to sort -tracer indices for more optimal access (eg contiguous indices), -prune the groups, and finalize group members. This should only be -called at the end of the `Tracers::init` so is defined as a -private function: ```c++ -private: - static int organize(); +static I4 getNumTracers(); ``` -An error code is returned if any error conditions are encountered. -#### 4.2.5 getNumTracers +#### 4.2.4 `getIndex` + +This function retrieves the tracer index based on the tracer name. -Retrieves the total number of selected tracers activated for the -current simulation. ```c++ -static int getNumTracers(); +// Retrieves tracer index from tracer name +static I4 getIndex(I4 &TracerIndex, ///< [out] Tracer index + const std::string &TracerName ///< [in] Tracer name + +); ``` -#### 4.2.6 getIndex +#### 4.2.5 `getAll` + +This function retrieves the full array of all tracers at a specified time +level. The current time level is specified by an integer value of 0, the +previous time level is -1, and so on. If the requested time level does not +exist, the function returns a negative integer. -This function retrieves the tracer index given a tracer name. ```c++ -static int getIndex(const std::string &TracerName); +// get a device array for all tracers +static I4 getAll(Array3DReal &TracerArray, ///< [out] tracer device array + const I4 TimeLevel ///< [in] time level index +); ``` -#### 4.2.7 getAll +#### 4.2.6 `getByIndex` (Single Tracer) + +This function retrieves a single tracer array by index and time level. +Similar to `getAll`, the current time level is specified by an integer value +of 0, and the previous time level is -1. If the requested time level does +not exist, the function returns a negative integer. -This function retrieves the full array of all tracers at a -given time level. ```c++ -static Array3DReal getAll(int TimeIndx ///< [in] time level index +// get a device array by tracer index +static I4 getByIndex(Array2DReal &TracerArray, ///< [out] tracer device array + const I4 TimeLevel, ///< [in] time level index + const I4 TracerIndex ///< [in] global tracer index ); ``` -#### 4.2.8 get (single tracer and time level) +#### 4.2.7 `getGroupRange` + +This function retrieves a pair of the group's start index and group length. -To retrieve a single tracer array by index and time level, we -provide the function: ```c++ -static Array2dReal get(int TrcrIndx, ///< [in] global tracer index - int TimeIndx ///< [in] time level index +// Retrieves a pair of (group start index, group length) +static I4 +getGroupRange(std::pair &GroupRange, ///< [out] Group range + const std::string &GroupName ///< [in] Group name + ); ``` -#### 4.2.9 getGroup +A typical use of this function might look like: -This function retrieves a vector of tracer indices that belong -to a group. The group name must be provided. ```c++ -static std::vector getGroup(const std::string &GroupName); +// Get group range +std::pair GroupRange; +I4 Err = Tracers::getGroupRange(GroupRange, GroupName); + +// Unpack group range +auto [StartIndex, GroupLength] = GroupRange; + +// Get all tracers at the current time level (0) +Array3DReal TracerArray; +Err = OMEGA::Tracers::getAll(TracerArray, 0); +if (Err != 0) + LOG_ERROR("getAll returns an error code: {}.", Err); + +OMEGA::parallelFor( + "ComputeGroupTendency", + {GroupLength, NCells, NVertLayers}, + KOKKOS_LAMBDA(Int iIndex, Int iCell, Int iVert) { + int iTracer = TracerArray[iIndex + StartIndex]; + // Perform operations on TracerArray(iTracer, iCell, iVert); + }); ``` -A typical use of this function might look like: -```c++ - Array4DReal TrcrArray = OMEGA::Tracers::getAll(); - std::vector EcoTracerIndices = - OMEGA::Tracers::getGroup("Ecosys") - - int NGrpMembers = EcoTracerIndices.size(); - OMEGA::parallelFor( - "ComputeGroupTendency", - {NGrpMembers, NCells, NVertLayers}, - KOKKOS_LAMBDA(Int iTracer, Int iCell, Int k) { - int TrcrIndx = EcoTracerIndices[iTracer]; - do stuff with TrcrArray(TrcrIndx, iCell, k); - }); -``` +#### 4.2.8 `isGroupMemberByIndex` -#### 4.2.10 isMember +When looping over tracers, this function returns `true` if the current +tracer index belongs to a group with a given name. -When looping over tracers, this function returns true if the -current tracer index belongs to a group of a given name. -```c++ -static bool isMember( - int Index, ///< [in] global tracer index - const std::string &GroupName); -``` -This provides a complement to the getGroup function so the -loop in the previous section could also look like: ```c++ - int NTracers = OMEGA::Tracers::getNumTracers; - Array4DReal TrcrArray = OMEGA::Tracers::getAll(); - - OMEGA::parallelFor( - "ComputeGroupTendency", - {NTracers, NCells, NVertLayers}, - KOKKOS_LAMBDA(Int iTracer, Int iCell, Int k) { - if (isMember(iTracer,"GrpName"){ - do stuff with TrcrArray(TrcrIndx, iCell, k); - } - }); +// Checks if a tracer is a member of a group by tracer index +static bool +isGroupMemberByIndex(const I4 TracerIndex, ///< [in] Tracer index + const std::string GroupName ///< [in] Group name +); ``` -#### 4.2.11 Managing time levels +#### 4.2.9 Managing Time Levels + +The exact process for managing time integration is not yet clear. We +assume a time integration module will assign and swap time indices (e.g., +CurTime, PrevTime) to point to the appropriate array index in all +prognostic arrays like tracers. The current time is represented by +an integer value of "0", while the previous time is represented as "-1", +and so on. -It is not clear yet how time integration will be managed. We assume -here that a time integration module will be assigning and swapping -time indices (eg CurTime, NewTime, OldTime?, etc.) to point to -appropriate array index in all prognostic arrays like tracers. -If these need to be managed within the tracer module, we may need -a function like swapTimeLevels to identify and modify time level -indices. +#### 4.2.10 `clear` -#### 4.2.12 clear +This function clears all defined tracers and cleans up memory. It returns +an error code if it is unable to deallocate or destroy all memory. -A function is needed to clear all defined tracers and clean up -memory. An error code is returned if it is unable to deallocate -or destroy all memory. ```c++ +// Clears all defined tracers and deallocates memory. static int clear(); ``` +Here is the corrected version with proper grammar: + +--- + ## 5 Verification and Testing ### 5.1 Test All -A sample TracerDefs.inc file will be created in the test directory -with the base tracers (T,S) and two tracers in a Debug group. -We will assume two time levels (Cur, New). A unit test driver -will read an input configuration that requests these two -tracer groups. The tracers will be initialized with artificial -data and basic arithmetic will be performed on the device. In one -loop nest, we will perform arithmetic on one group using the -vector of group tracer indices. In a second loop nest, we will -perform operations on the second group using the isMembers function -to select. - - tests requirements 2.1-2.4, 2.6 - -All other requirements will need to be tested in later system -tests once IOStreams and time integration schemes have been -implemented. +A sample `TracerDefs.inc` file will be located in the source directory, +containing the base tracers and three tracers in a Debug group. We will +assume two time levels (Cur: 0, Prev: -1). A unit test driver will read an +input configuration that requests these two tracer groups. The tracers will +be initialized with artificial data, and basic arithmetic will be performed +on the device. In one loop, we will perform arithmetic on one group using +the vector of group tracer indices. In a second loop, we will perform +operations on the second group using the `isMembers` function to select tracers. + +This test covers requirements 2.1–2.4 and 2.6. + +All other requirements will need to be tested in later system tests once +I/O streams and time integration schemes have been implemented. diff --git a/components/omega/doc/devGuide/Tracers.md b/components/omega/doc/devGuide/Tracers.md new file mode 100644 index 000000000000..384c00c24bd0 --- /dev/null +++ b/components/omega/doc/devGuide/Tracers.md @@ -0,0 +1,210 @@ +(omega-dev-tracers)= + +# Tracers + +Tracers refer to either heat (temperature) or material carried by a fluid +parcel (e.g., salt, chemical, or biological constituents). + +To manage tracer definitions, users will update the `TracerDefs.inc` file +located in the `omega/src/ocn` directory of E3SM. The file contains +`define` C++ function calls, each of which defines a tracer as shown below: + +```c++ +define( + "Temp", // Name of variable + "Potential Temperature", // Long name or description + "degree_C", // Units + "sea_water_potential_temperature", // CF standard name + -273.15, // Min valid field value + 100.0, // Max valid field value + 1.e33 // Fill value for undefined entries +); +``` + +To add a new tracer, simply call the `define` function with the appropriate +arguments. + +The following sections explain the concept and implementation of OMEGA tracers. + +## Concepts + +### Arrays in Host and Device Memory + +All tracers are stored in a vector of 3-dimensional device and host arrays. + +```c++ + std::vector TracerArrays; // Device: multiple time levels + std::vector TracerArraysH; // Host: multiple time levels +``` + +Each device and host array in the vector has dimensions corresponding to the +number of tracers, the number of cells in the rank, and the vertical levels. + +The host tracer array (`TracerArraysH`) is internally managed, so users +should not directly modify it in most cases. + +### Tracer Groups + +Each tracer is assigned to a tracer group defined in the OMEGA YAML +configuration file, such as `Default.yml`. A tracer should be assigned to only one group. + +To access the member tracers of a group in the code, users can use the +`getGroupRange` function as shown below: + +```c++ +// Retrieves a pair of (group start index, group length) +static I4 +getGroupRange(std::pair &GroupRange, ///< [out] Group range + const std::string &GroupName ///< [in] Group name +); +``` + +A typical use of this function might look like: + +```c++ +// Get group range +std::pair GroupRange; +I4 Err = Tracers::getGroupRange(GroupRange, GroupName); + +// Unpack group range +auto [StartIndex, GroupLength] = GroupRange; + +// Get all tracers at the current time level (0) +Array3DReal TracerArray; +Err = OMEGA::Tracers::getAll(TracerArray, 0); +if (Err != 0) + LOG_ERROR("getAll returns an error code: {}.", Err); + + +OMEGA::parallelFor( + "ComputeGroupTendency", + {GroupLength, NCells, NVertLayers}, + KOKKOS_LAMBDA(Int iIndex, Int iCell, Int iVert) { + int iTracer = TracerArray[iIndex + StartIndex]; + // Perform operations on TracerArray(iTracer, iCell, iVert); + }); +``` + +### Time Levels + +During the initialization of Tracers, the number of time levels is determined, +and the length of the tracer vectors (`TracerArrays` and `TracerArraysH`) +is set accordingly. The Tracers class internally manages the index of the +current time level in the arrays. To access a specific time level in the +arrays, users will use the integer "0" or negative integers. For example, +the following code retrieves the device tracer arrays for the current and +previous time levels: + +```c++ +// Get all tracers at the current time level (0) +Array3DReal CurrentTracerArray; +I4 Err1 = OMEGA::Tracers::getAll(CurrentTracerArray, 0); + +// Get all tracers at the previous time level (-1) +Array3DReal PreviousTracerArray; +I4 Err2 = OMEGA::Tracers::getAll(PreviousTracerArray, -1); +``` + +## Initialization and Finalization + +To use Tracers, the class should first be initialized by calling the +`init()` function and finalized by calling `finalize()`. These function +calls are typically handled within the initialization and finalization of +OMEGA itself, so users may not need to call them separately. + +## Key APIs + +### `getAll` and `getAllHost` + +These functions return all device and host tracer arrays, respectively. If +the specified `TimeLevel` does not exist, they return a negative integer. + +```c++ +static HostArray3DReal getAllHost( + HostArray3DReal &TracerArrayH, ///< [out] tracer host array + const I4 TimeLevel ///< [in] Time level index +); + +static Array3DReal getAll( + Array3DReal &TracerArray, ///< [out] tracer device array + const I4 TimeLevel ///< [in] Time level index +); +``` + +### `getGroupRange` + +`getGroupRange` returns a pair of two `I4` values representing the group +start index and the group length. If the group is not found, it returns +negative integers. + +```c++ +static I4 getGroupRange( + std::pair &GroupRange, ///< [out] Group range + const std::string &GroupName ///< [in] Group name +); +``` + +### `getByIndex` and `getHostByIndex` + +These functions return device and host tracer arrays, respectively, based +on the `TracerIndex`. If the specified `TimeLevel` and/or `TracerIndex` +does not exist, they return a negative integer. + +```c++ +static Array2DReal getByIndex( + Array2DReal &TracerArray, ///< [out] tracer device array + const I4 TimeLevel, ///< [in] Time level index + const I4 TracerIndex ///< [in] Global tracer index +); + +static HostArray2DReal getHostByIndex( + HostArray2DReal &TracerArrayH, ///< [out] tracer host array + const I4 TimeLevel, ///< [in] Time level index + const I4 TracerIndex ///< [in] Global tracer index +); +``` + +### `getIndex` + +`getIndex` returns the index of the tracer specified by the `TracerName` +argument. If the tracer is not found, it returns a negative integer. + +```c++ +static I4 getIndex( + I4 &TracerIndex, ///< [out] Tracer index + const std::string &TracerName ///< [in] Tracer name +); +``` + +### `getFieldByIndex` + +`getFieldByIndex` returns the `Field` object associated with the tracer +specified by the `TracerIndex`. If not found, it returns `nullptr`. + +```c++ +// Returns a field by tracer index. If it does not exist, returns nullptr +static std::shared_ptr +getFieldByIndex(const I4 TracerIndex ///< [in] Global tracer index +); +``` + +### `updateTimeLevels` + +`updateTimeLevels` increments the current time level by one. If the current +time level exceeds the number of time levels, it returns to `0`. It also +exchanges the halo cells of all tracers at the current time level and updates +the associated field with the new tracer arrays. + +```c++ +/// Increment time levels +static I4 updateTimeLevels(); +``` + +### `getNumTracers` + +`getNumTracers` returns the total number of tracers used in the simulation. + +```c++ +// Get total number of tracers +static I4 getNumTracers(); +``` diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 7d4f5cf507e3..240c51e70f69 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -43,6 +43,7 @@ userGuide/OceanState userGuide/TimeMgr userGuide/TimeStepping userGuide/Reductions +userGuide/Tracers ``` ```{toctree} @@ -76,6 +77,7 @@ devGuide/OceanState devGuide/TimeMgr devGuide/TimeStepping devGuide/Reductions +devGuide/Tracers ``` ```{toctree} diff --git a/components/omega/doc/userGuide/Tracers.md b/components/omega/doc/userGuide/Tracers.md new file mode 100644 index 000000000000..09c83983c2e2 --- /dev/null +++ b/components/omega/doc/userGuide/Tracers.md @@ -0,0 +1,42 @@ +(omega-user-tracers)= + +# Tracers + +Tracers refer to either heat (temperature) or material carried by a fluid +parcel (e.g., salt, chemical, or biological constituents). + +To manage tracer definitions, users will update the `TracerDefs.inc` file +located in the `omega/src/ocn` directory of E3SM. The file contains +`define` C++ function calls, each of which defines a tracer as shown below: + +```c++ +define( + "Temp", // Name of variable + "Potential Temperature", // Long name or description + "degree_C", // Units + "sea_water_potential_temperature", // CF standard name + -273.15, // Min valid field value + 100.0, // Max valid field value + 1.e33 // Fill value for undefined entries +); +``` + +To add a new tracer, simply call the `define` function with the appropriate +arguments. + +Note that not all tracers defined in `TracerDefs.inc` will be used during +a simulation. To select tracers and groups of tracers, users will configure +YAML files in OMEGA, such as `Default.yml` in the `omega/configs` directory, +as shown below: + +```yaml +omega: + Tracers: + Base: [Temp, Salt] + Debug: [Debug1, Debug2, Debug3] + [other individual tracers or groups as needed] +``` + +In the above example, two tracer groups (Base and Debug) are selected. The +Base group includes the `Temp` and `Salt` tracers, while the Debug group +includes `Debug1`, `Debug2`, and `Debug3`. diff --git a/components/omega/src/ocn/TracerDefs.inc b/components/omega/src/ocn/TracerDefs.inc new file mode 100644 index 000000000000..278a3823dc86 --- /dev/null +++ b/components/omega/src/ocn/TracerDefs.inc @@ -0,0 +1,25 @@ +//===-- ocn/TracerDefs.inc --------------------*- C++ -*-===// +// +/// \file +/// \brief OMEGA Tracer Definition File +/// +/// This file defines all OMEGA tracers. It is imported using +/// the C++ #include preprocessor directive in Tracers.cpp. Only the tracer +/// definitions selected in the Tracers section of the OMEGA YAML configuration +/// file will be allocated into memory. +//===----------------------------------------------------------------------===// + +define("Temp", ///< Name of tracer + "Potential Temperature", ///< Long name or description + "degree_C", ///< Units + "sea_water_potential_temperature", ///< CF standard Name + -273.15, ///< min valid field value + 100.0, ///< max valid field value + 1.e33 ///< value for undef entries +); + +define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0, 1.e33); + +define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, 1.e33); +define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, 1.e33); +define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, 1.e33); diff --git a/components/omega/src/ocn/Tracers.cpp b/components/omega/src/ocn/Tracers.cpp new file mode 100644 index 000000000000..d627a49b15c9 --- /dev/null +++ b/components/omega/src/ocn/Tracers.cpp @@ -0,0 +1,527 @@ +//===-- OMEGA Tracers Implementation -----------------------------*- C++ +//-*-===// +// +/// \file +/// \brief OMEGA Tracers Implementation +/// +// +//===-----------------------------------------------------------------------===/ + +#include "Tracers.h" +#include "Decomp.h" +#include "IO.h" +#include "Logging.h" +#include "TimeStepper.h" + +#include + +namespace OMEGA { + +// Initialize static member variables +std::vector Tracers::TracerArrays; +std::vector Tracers::TracerArraysH; + +std::map> Tracers::TracerGroups; +std::map Tracers::TracerIndexes; +std::map Tracers::TracerNames; +std::vector Tracers::TracerDimNames = {"NCells", "NVertLevels"}; + +Halo *Tracers::MeshHalo = nullptr; + +I4 Tracers::NCellsOwned = 0; +I4 Tracers::NCellsAll = 0; +I4 Tracers::NCellsSize = 0; +I4 Tracers::NTimeLevels = 0; +I4 Tracers::NVertLevels = 0; +I4 Tracers::CurTimeIndex = 0; +I4 Tracers::NumTracers = 0; + +//--------------------------------------------------------------------------- +// Internal Utilities +//--------------------------------------------------------------------------- +std::string Tracers::packTracerFieldName(const std::string &TracerName) { + return "Tracer" + TracerName; +} + +//--------------------------------------------------------------------------- +// Initialization +//--------------------------------------------------------------------------- +I4 Tracers::init() { + + I4 Err = 0; + + LOG_INFO("Tracers::init() called"); + + // Retrieve mesh cell/edge/vertex totals from Decomp + HorzMesh *DefHorzMesh = HorzMesh::getDefault(); + + NCellsOwned = DefHorzMesh->NCellsOwned; + NCellsAll = DefHorzMesh->NCellsAll; + NCellsSize = DefHorzMesh->NCellsSize; + NVertLevels = DefHorzMesh->NVertLevels; + + MeshHalo = Halo::getDefault(); + + auto *DefTimeStepper = TimeStepper::getDefault(); + if (!DefTimeStepper) { + LOG_ERROR("TimeStepper needs to be initialized before Tracers"); + return -1; + } + NTimeLevels = DefTimeStepper->getNTimeLevels(); + + if (NTimeLevels < 2) { + LOG_ERROR("Tracers: the number of time level is lower than 2"); + return -2; + } + + CurTimeIndex = 0; + + // load Tracers configs + Config *OmegaConfig = Config::getOmegaConfig(); + Config TracersConfig("Tracers"); + Err = OmegaConfig->get(TracersConfig); + if (Err != 0) { + LOG_ERROR("Tracers: Tracers group not found in Config"); + return -3; + } + + NumTracers = 0; + I4 TracerIndex = 0; + + // get tracers group and tracer names + for (auto It = TracersConfig.begin(); It != TracersConfig.end(); ++It) { + + I4 GroupStartIndex = TracerIndex; + + std::string GroupName; + I4 GroupNameErr = OMEGA::Config::getName(It, GroupName); + if (GroupNameErr != 0) { + LOG_ERROR("Tracers: {} tracer group name not found in TracersConfig", + GroupName); + return -4; + } + + std::vector _TracerNames; + I4 TracerNamesErr = TracersConfig.get(GroupName, _TracerNames); + if (TracerNamesErr != 0) { + LOG_ERROR("Tracers: {} group tracers not found in TracersConfig", + GroupName); + return -5; + } + + for (auto _TracerName : _TracerNames) { + TracerIndexes[_TracerName] = TracerIndex; + TracerIndex++; + } + + TracerGroups[GroupName] = + std::pair(GroupStartIndex, TracerIndex - GroupStartIndex); + FieldGroup::create("TracerGroup" + GroupName); + } + + // total number of tracers + NumTracers = TracerIndex; + + // Initialize tracers arrays for device and host + TracerArrays.resize(NTimeLevels); + TracerArraysH.resize(NTimeLevels); + + // Allocate tracers data array and assign to tracers arrays + for (I4 TimeIndex = 0; TimeIndex < NTimeLevels; ++TimeIndex) { + TracerArrays[TimeIndex] = + Array3DReal("TracerTime" + std::to_string(TimeIndex), NumTracers, + NCellsSize, NVertLevels); + TracerArraysH[TimeIndex] = + HostArray3DReal("TracerHTime" + std::to_string(TimeIndex), NumTracers, + NCellsSize, NVertLevels); + } + +// Read tracer definitions from file +#include "TracerDefs.inc" + + // Check if all tracers defined in config file are loaded + if (TracerIndexes.size() != TracerNames.size()) { + LOG_ERROR("Tracer: not all tracers defined in config file is loaded."); + return -6; + } + + // Add Fields to FieldGroup + for (const auto &GroupPair : TracerGroups) { + auto GroupName = GroupPair.first; + std::string TracerFieldGroupName = "TracerGroup" + GroupName; + auto TracerFieldGroup = FieldGroup::get(TracerFieldGroupName); + + std::vector _TracerNames; + TracersConfig.get(GroupName, _TracerNames); + + for (auto _TracerName : _TracerNames) { + std::string TracerFieldName = packTracerFieldName(_TracerName); + + // add tracer Field to field group + Err = TracerFieldGroup->addField(TracerFieldName); + if (Err != 0) { + LOG_ERROR("Error adding {} to field group {}", TracerFieldName, + TracerFieldGroupName); + return -7; + } + + // Associate Field with data + I4 TracerIndex = TracerIndexes[_TracerName]; + std::shared_ptr TracerField = Field::get(TracerFieldName); + + // Create a 2D subview by fixing the first dimension (TracerIndex) + Array2DReal TracerSubview = Kokkos::subview( + TracerArrays[CurTimeIndex], TracerIndex, Kokkos::ALL, Kokkos::ALL); + Err = TracerField->attachData(TracerSubview); + if (Err != 0) { + LOG_ERROR("Error attaching data array to field {}", + TracerFieldName); + return -8; + } + } + } + + return 0; +} + +//--------------------------------------------------------------------------- +// Define tracers +//--------------------------------------------------------------------------- +I4 Tracers::define(const std::string &Name, const std::string &Description, + const std::string &Units, const std::string &StdName, + const Real ValidMin, const Real ValidMax, + const Real FillValue) { + + // Do nothing if this tracer is not selected + if (TracerIndexes.find(Name) == TracerIndexes.end()) { + return 0; + } + + auto TracerIndex = TracerIndexes[Name]; + + // Return error if tracer already exists + if (TracerNames.find(TracerIndex) != TracerNames.end()) { + LOG_ERROR("Tracers: Tracer '{}' already exists", Name); + return -1; + } + + // set tracer index to name mapping + TracerNames[TracerIndex] = Name; + + // create a tracer field + std::string TracerFieldName = packTracerFieldName(Name); + auto TracerField = Field::create(TracerFieldName, Description, Units, + StdName, ValidMin, ValidMax, FillValue, + TracerDimNames.size(), TracerDimNames); + if (!TracerField) { + LOG_ERROR("Tracers: Tracer field '{}' is not created", TracerFieldName); + return -2; + } + + return 0; +} + +//--------------------------------------------------------------------------- +// Deallocate tracer arrays +//--------------------------------------------------------------------------- +I4 Tracers::clear() { + + LOG_INFO("Tracers::clear() called"); + + // Deallocate memory for tracer arrays + TracerArrays.clear(); + TracerArraysH.clear(); + + TracerGroups.clear(); + TracerIndexes.clear(); + TracerNames.clear(); + + NumTracers = 0; + NCellsOwned = 0; + NCellsAll = 0; + NCellsSize = 0; + NTimeLevels = 0; + NVertLevels = 0; + + return 0; +} + +//--------------------------------------------------------------------------- +// Query tracers +//--------------------------------------------------------------------------- + +I4 Tracers::getNumTracers() { return NumTracers; } + +I4 Tracers::getIndex(I4 &TracerIndex, const std::string &TracerName) { + // Check if tracer exists + if (TracerIndexes.find(TracerName) != TracerIndexes.end()) { + TracerIndex = TracerIndexes[TracerName]; + return 0; // Success + } + + LOG_ERROR("Tracers: Tracer index for '{}' is not found.", TracerName); + return -1; // Tracer not found +} + +I4 Tracers::getName(std::string &TracerName, const I4 TracerIndex) { + if (TracerNames.find(TracerIndex) != TracerNames.end()) { + TracerName = TracerNames[TracerIndex]; + return 0; // Success + } + + LOG_ERROR("Tracers: Tracer name for index '{}' is not found.", TracerIndex); + return -1; // Tracer index not found +} + +//--------------------------------------------------------------------------- +// get Tracer arrays +// TimeLevel == [0:current, -1:previous, -2:two times ago, ...] +//--------------------------------------------------------------------------- +I4 Tracers::getAll(Array3DReal &TracerArray, const I4 TimeLevel) { + I4 Err = 0; + + if (TimeLevel <= 0 && (TimeLevel + NTimeLevels) > 0) { + I4 TimeIndex = (TimeLevel + CurTimeIndex + NTimeLevels) % NTimeLevels; + TracerArray = TracerArrays[TimeIndex]; + } else { + LOG_ERROR("Tracers: Time level {} is out of range", TimeLevel); + Err = -1; + } + + return Err; +} + +I4 Tracers::getByIndex(Array2DReal &TracerArray, const I4 TimeLevel, + const I4 TracerIndex) { + // Check if time level is valid + if (TimeLevel > 0 || (TimeLevel + NTimeLevels) <= 0) { + LOG_ERROR("Tracers: Time level {} is out of range", TimeLevel); + return -1; + } + + // Check if tracer index is valid + if (TracerIndex < 0 || TracerIndex >= NumTracers) { + LOG_ERROR("Tracers: Tracer index {} is out of range", TracerIndex); + return -2; + } + + I4 TimeIndex = (TimeLevel + CurTimeIndex + NTimeLevels) % NTimeLevels; + TracerArray = Kokkos::subview(TracerArrays[TimeIndex], TracerIndex, + Kokkos::ALL, Kokkos::ALL); + return 0; +} + +I4 Tracers::getByName(Array2DReal &TracerArray, const I4 TimeLevel, + const std::string &TracerName) { + // Check if tracer exists + if (TracerIndexes.find(TracerName) == TracerIndexes.end()) { + LOG_ERROR("Tracers: Tracer '{}' does not exist", TracerName); + return -1; + } + + int Err = getByIndex(TracerArray, TimeLevel, TracerIndexes[TracerName]); + if (Err != 0) + return -2; + + return 0; +} + +I4 Tracers::getAllHost(HostArray3DReal &TracerArrayH, const I4 TimeLevel) { + + I4 Err = 0; + + if (TimeLevel <= 0 && (TimeLevel + NTimeLevels) > 0) { + I4 TimeIndex = (TimeLevel + CurTimeIndex + NTimeLevels) % NTimeLevels; + TracerArrayH = TracerArraysH[TimeIndex]; + } else { + LOG_ERROR("Tracers: Time level {} is out of range", TimeLevel); + Err = -1; + } + return Err; +} + +I4 Tracers::getHostByIndex(HostArray2DReal &TracerArrayH, const I4 TimeLevel, + const I4 TracerIndex) { + // Check if time level is valid + if (TimeLevel > 0 || (TimeLevel + NTimeLevels) <= 0) { + LOG_ERROR("Tracers: Time level {} is out of range", TimeLevel); + return -1; + } + + // Check if tracer index is valid + if (TracerIndex < 0 || TracerIndex >= NumTracers) { + LOG_ERROR("Tracers: Tracer index {} is out of range", TracerIndex); + return -2; + } + + I4 TimeIndex = (TimeLevel + CurTimeIndex + NTimeLevels) % NTimeLevels; + TracerArrayH = Kokkos::subview(TracerArraysH[TimeIndex], TracerIndex, + Kokkos::ALL, Kokkos::ALL); + return 0; +} + +I4 Tracers::getHostByName(HostArray2DReal &TracerArrayH, const I4 TimeLevel, + const std::string &TracerName) { + // Check if tracer exists + if (TracerIndexes.find(TracerName) == TracerIndexes.end()) { + LOG_ERROR("Tracers: Tracer '{}' does not exist", TracerName); + return -1; + } + + // Get the index of the tracer + I4 TracerIndex = TracerIndexes[TracerName]; + I4 Err = getHostByIndex(TracerArrayH, TimeLevel, TracerIndex); + if (Err != 0) + return -2; + + return 0; +} + +//--------------------------------------------------------------------------- +// get Fields +//--------------------------------------------------------------------------- +std::shared_ptr Tracers::getFieldByName(const std::string &TracerName) { + // Check if tracer exists + if (TracerIndexes.find(TracerName) == TracerIndexes.end()) { + LOG_ERROR("Tracers: Tracer '{}' does not exist", TracerName); + return nullptr; + } + + return Field::get(packTracerFieldName(TracerName)); +} + +std::shared_ptr Tracers::getFieldByIndex(const I4 TracerIndex) { + // Check if tracer index is valid + if (TracerIndex < 0 || TracerIndex >= NumTracers) { + LOG_ERROR("Tracers: Tracer index {} is out of range", TracerIndex); + return nullptr; + } + + return getFieldByName(TracerNames[TracerIndex]); +} + +//--------------------------------------------------------------------------- +// Query tracer group +//--------------------------------------------------------------------------- +std::vector Tracers::getGroupNames() { + std::vector GroupNames; + + for (const auto &GroupPair : TracerGroups) { + GroupNames.push_back(GroupPair.first); + } + + return GroupNames; +} + +I4 Tracers::getGroupRange(std::pair &GroupRange, + const std::string &GroupName) { + auto it = TracerGroups.find(GroupName); + if (it != TracerGroups.end()) { + GroupRange = it->second; + return 0; + } + + return -1; +} + +bool Tracers::isGroupMemberByIndex(const I4 TracerIndex, + const std::string GroupName) { + auto it = TracerGroups.find(GroupName); + if (it != TracerGroups.end()) { + I4 StartIndex = it->second.first; + I4 GroupLength = it->second.second; + return TracerIndex >= StartIndex && + TracerIndex < StartIndex + GroupLength; + } + + return false; +} + +bool Tracers::isGroupMemberByName(const std::string &TracerName, + const std::string &GroupName) { + + I4 TracerIndex; + if (getIndex(TracerIndex, TracerName) != 0) { + return false; + } + + return isGroupMemberByIndex(TracerIndex, GroupName); +} + +//--------------------------------------------------------------------------- +// deep copy at a time level +// TimeLevel == [0:current, -1:previous, -2:two times ago, ...] +//--------------------------------------------------------------------------- +I4 Tracers::copyToDevice(const I4 TimeLevel) { + + // Check if time level is valid + if (TimeLevel > 0 || (TimeLevel + NTimeLevels) <= 0) { + LOG_ERROR("Tracers: Time level {} is out of range", TimeLevel); + return -1; + } + + I4 TimeIndex = (TimeLevel + CurTimeIndex + NTimeLevels) % NTimeLevels; + deepCopy(TracerArrays[TimeIndex], TracerArraysH[TimeIndex]); + + return 0; +} + +I4 Tracers::copyToHost(const I4 TimeLevel) { + + // Check if time level is valid + if (TimeLevel > 0 || (TimeLevel + NTimeLevels) <= 0) { + LOG_ERROR("Tracers: Time level {} is out of range", TimeLevel); + return -1; + } + + I4 TimeIndex = (TimeLevel + CurTimeIndex + NTimeLevels) % NTimeLevels; + deepCopy(TracerArraysH[TimeIndex], TracerArrays[TimeIndex]); + + return 0; +} + +//--------------------------------------------------------------------------- +// halo exchange +// TimeLevel == [0:current, -1:previous, -2:two times ago, ...] +//--------------------------------------------------------------------------- +I4 Tracers::exchangeHalo(const I4 TimeLevel) { + // TODO: copy only halo cells + copyToHost(TimeLevel); + + I4 TimeIndex = (TimeLevel + CurTimeIndex + NTimeLevels) % NTimeLevels; + MeshHalo->exchangeFullArrayHalo(TracerArraysH[TimeIndex], OnCell); + + // TODO: copy only halo cells + copyToDevice(TimeLevel); +} + +//--------------------------------------------------------------------------- +// update time level +//--------------------------------------------------------------------------- +I4 Tracers::updateTimeLevels() { + + // Exchange halo + exchangeHalo(0); + + CurTimeIndex = (CurTimeIndex + 1) % NTimeLevels; + + // Update TracerField data associations + for (const auto &TracerPair : TracerIndexes) { + auto TracerFieldName = packTracerFieldName(TracerPair.first); + auto TracerIndex = TracerPair.second; + + std::shared_ptr TracerField = Field::get(TracerFieldName); + + Array2DReal TracerSubview = Kokkos::subview( + TracerArrays[CurTimeIndex], TracerIndex, Kokkos::ALL, Kokkos::ALL); + I4 Err = TracerField->attachData(TracerSubview); + if (Err != 0) { + LOG_ERROR("Error attaching data array to field {}", TracerFieldName); + return Err; + } + } + + return 0; +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/Tracers.h b/components/omega/src/ocn/Tracers.h new file mode 100644 index 000000000000..31b77fc8b72f --- /dev/null +++ b/components/omega/src/ocn/Tracers.h @@ -0,0 +1,212 @@ +#ifndef OMEGA_TRACERS_H +#define OMEGA_TRACERS_H +//===-- ocn/Tracers.h - tracers --------------------*- C++ -*-===// +// +/// \file +/// \brief Define Tracers class +/// +/// Tracers class define tracer arrays for host and device memories and +/// methods to handle the arrays. All of variables and methods are static +/// because once tracers are created, they never change. +//===----------------------------------------------------------------------===// + +#include "DataTypes.h" +#include "Decomp.h" +#include "Field.h" +#include "Halo.h" + +namespace OMEGA { + +/// The Tracers class provides a container for tracer arrays and methods +class Tracers { + + private: + // total number of tracers defined at intialization + static I4 NumTracers; + + // static storage of the tracer arrays for device and host memory spaces + // The arrays are 3 dimensional arrays with Tracer, Cell, and Vertical + // dimensions in order + static std::vector + TracerArrays; ///< TimeLevels -> [Tracer, Cell, Vert] + static std::vector + TracerArraysH; ///< TimeLevels -> [Tracer, Cell, Vert] + + // maps for managing tracer groups + // Key of this map is a group name and + // Value is a pair of GroupStartIndex and GroupLength + static std::map> TracerGroups; + + // maps for matching tracer names with indices (both directions) + static std::map TracerIndexes; + static std::map TracerNames; + + // for halo exchange + static Halo *MeshHalo; + + // Tracer dimension names + static std::vector TracerDimNames; + + // Current time index + // this index is circular so that it returns to index 0 + // if it is over max index + static I4 CurTimeIndex; ///< Time dimension array index for current level + + // pack tracer field name + static std::string packTracerFieldName(const std::string &TracerName); + + // locally defines all tracers but do not allocates memory + static I4 + define(const std::string &Name, ///< [in] Name of tracer + const std::string &Description, ///< [in] Long name or description + const std::string &Units, ///< [in] Units + const std::string &StdName, ///< [in] CF standard Name + const Real ValidMin, ///< [in] min valid field value + const Real ValidMax, ///< [in] max valid field value + const Real FillValue ///< [in] value for undef entries + ); + + public: + static I4 NTimeLevels; ///< Number of time levels in tracer variable arrays + static I4 + NVertLevels; ///< Number of vertical levels in tracer variable arrays + static I4 NCellsOwned; ///< Number of cells owned by this task + static I4 NCellsAll; ///< Total number of local cells (owned + all halo) + static I4 + NCellsSize; ///< Array size (incl padding, bndy cell) for cell arrays + + //--------------------------------------------------------------------------- + // Initialization + //--------------------------------------------------------------------------- + /// read tracer defintions, allocate tracer arrays and initializes the + /// tracers + static I4 init(); + + /// deallocates tracer arrays + static I4 clear(); + + //--------------------------------------------------------------------------- + // Query tracers + //--------------------------------------------------------------------------- + + // get total number of tracers + static I4 getNumTracers(); + + // get tracer index from tracer name + static I4 getIndex(I4 &TracerIndex, ///< [out] tracer index + const std::string &TracerName ///< [in] tracer name + ); + + // get tracer name from tracer index + static I4 getName(std::string &TracerName, ///< [out] tracer name + const I4 TracerIndex ///< [in] tracer index + ); + + // get a device array for all tracers + static I4 getAll(Array3DReal &TracerArray, ///< [out] tracer device array + const I4 TimeLevel ///< [in] time level index + ); + + // get a device array by tracer index + static I4 getByIndex(Array2DReal &TracerArray, ///< [out] tracer device array + const I4 TimeLevel, ///< [in] time level index + const I4 TracerIndex ///< [in] global tracer index + ); + + // get a device array by tracer name + static I4 + getByName(Array2DReal &TracerArray, ///< [out] tracer device array + const I4 TimeLevel, ///< [in] time level index + const std::string &TracerName ///< [in] global tracer name + ); + + // get a host array for all tracers + static I4 + getAllHost(HostArray3DReal &TracerArrayH, ///< [out] tracer host array + const I4 TimeLevel ///< [in] time level index + ); + + // get a host array by tracer index + static I4 + getHostByIndex(HostArray2DReal &TracerArrayH, ///< [out] tracer host array + const I4 TimeLevel, ///< [in] time level index + const I4 TracerIndex ///< [in] global tracer index + ); + + // get a host array by tracer name + static I4 + getHostByName(HostArray2DReal &TracerArrayH, ///< [out] tracer host array + const I4 TimeLevel, ///< [in] time level index + const std::string &TracerName ///< [in] global tracer name + ); + + // get a field by tracer index. If not found, returns nullptr + static std::shared_ptr + getFieldByIndex(const I4 TracerIndex ///< [in] global tracer index + ); + + // get a field by tracer name. If not found, returns nullptr + static std::shared_ptr + getFieldByName(const std::string &TracerName ///< [in] global tracer name + ); + + //--------------------------------------------------------------------------- + // Tracer group query + //--------------------------------------------------------------------------- + + // return a vector of group names. + static std::vector getGroupNames(); + + // get a pair of (group start index, group length) + static I4 getGroupRange(std::pair &GroupRange, ///< [out] group range + const std::string &GroupName ///< [in] group name + ); + + // check if a tracer is a member of group by tracer index + static bool + isGroupMemberByIndex(const I4 TracerIndex, ///< [in] tracer index + const std::string GroupName ///< [in] group name + ); + + // check if a tracer is a member of group by tracer name + static bool isGroupMemberByName( + const std::string &TracerName, ///< [in] global tracer name + const std::string &GroupName ///< [in] group name + ); + + //--------------------------------------------------------------------------- + // Halo exchange and update time level + //--------------------------------------------------------------------------- + + /// Exchange halo + static I4 exchangeHalo(const I4 TimeLevel ///< [in] tracer time level + ); + + /// increment time levels + static I4 updateTimeLevels(); + + //--------------------------------------------------------------------------- + // Device-Host data movement + //--------------------------------------------------------------------------- + + /// Copy tracers variables from host to device + static I4 copyToDevice(const I4 TimeLevel ///< [in] tracer time level + ); + + /// Copy tracers variables from device to host + static I4 copyToHost(const I4 TimeLevel ///< [in] tracer time level + ); + + //--------------------------------------------------------------------------- + // Forbid copy and move construction + //--------------------------------------------------------------------------- + + Tracers(const Tracers &) = delete; + Tracers(Tracers &&) = delete; + +}; // end class Tracers + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// +#endif // defined OMEGA_TRACERS_H diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 00606d486c17..ff9b5bbedbdd 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -338,3 +338,14 @@ add_omega_test( drivers/StandaloneDriverTest.cpp "-n;8" ) + +################## +# Tracers test +################## + +add_omega_test( + TRACERS_TEST + testTracers.exe + ocn/TracersTest.cpp + "-n;8" +) diff --git a/components/omega/test/ocn/TracersTest.cpp b/components/omega/test/ocn/TracersTest.cpp new file mode 100644 index 000000000000..03025cf42f2a --- /dev/null +++ b/components/omega/test/ocn/TracersTest.cpp @@ -0,0 +1,491 @@ +//===-- Test driver for OMEGA Tracers -----------------------------*- C++ +//-*-===/ +// +/// \file +/// \brief Test driver for OMEGA tracers class +/// +/// This driver tests that the OMEGA tracers module. +// +//===-----------------------------------------------------------------------===/ + +#include "Tracers.h" + +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "Field.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "IO.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OmegaKokkos.h" +#include "TimeStepper.h" +#include "mpi.h" + +#include + +using namespace OMEGA; + +// misc. variables for testing +const Real RefReal = 3.0; + +//------------------------------------------------------------------------------ +// The initialization routine for Tracers testing. It calls various +// init routines, including the creation of the default decomposition. + +I4 initTracersTest() { + + I4 Err = 0; + + // Initialize the Machine Environment class - this also creates + // the default MachEnv. Then retrieve the default environment and + // some needed data members. + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + initLogging(DefEnv); + + // Open config file + Config("Omega"); + Err = Config::readAll("omega.yml"); + if (Err != 0) { + LOG_ERROR("Tracers: Error reading config file"); + return Err; + } + + // Initialize the IO system + Err = IO::init(DefComm); + if (Err != 0) { + LOG_ERROR("Tracers: error initializing parallel IO"); + return Err; + } + + // Create the default decomposition (initializes the decomposition) + Err = Decomp::init(); + if (Err != 0) { + LOG_ERROR("Tracers: error initializing default decomposition"); + return Err; + } + + // Initialize the default halo + Err = Halo::init(); + if (Err != 0) { + LOG_ERROR("Tracers: error initializing default halo"); + return Err; + } + + // Initialize the default mesh + Err = HorzMesh::init(); + if (Err != 0) { + LOG_ERROR("Tracers: error initializing default mesh"); + return Err; + } + + // Initialize the default time stepper + Err = TimeStepper::init(); + if (Err != 0) { + LOG_ERROR("Tracers: error initializing default time stepper"); + return Err; + } + return 0; +} + +//------------------------------------------------------------------------------ +// The test driver for Tracers infrastructure +// +int main(int argc, char *argv[]) { + + int RetVal = 0; + I4 Err; + I4 count = 0; + + // Initialize the global MPI environment + MPI_Init(&argc, &argv); + Kokkos::initialize(); + { + + // Call initialization routine + Err = initTracersTest(); + if (Err != 0) + LOG_ERROR("Tracers: Error initializing"); + + // Get MPI vars if needed + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm Comm = DefEnv->getComm(); + I4 MyTask = DefEnv->getMyTask(); + I4 NumTasks = DefEnv->getNumTasks(); + bool IsMaster = DefEnv->isMasterTask(); + + HorzMesh *DefHorzMesh = HorzMesh::getDefault(); + Decomp *DefDecomp = Decomp::getDefault(); + Halo *DefHalo = Halo::getDefault(); + + // initialize Tracers infrastructure + Err = Tracers::init(); + if (Err != 0) { + RetVal += 1; + LOG_ERROR("Tracers: initialzation FAIL"); + } + + // Get group names + std::vector GroupNames = Tracers::getGroupNames(); + + // Check if "Base" group exists + if (std::find(GroupNames.begin(), GroupNames.end(), "Base") != + GroupNames.end()) { + LOG_INFO("Tracers: Group, 'Base', exists PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: Group, 'Base', does not exist FAIL"); + } + + // Check if "Debug" group exists + if (std::find(GroupNames.begin(), GroupNames.end(), "Debug") != + GroupNames.end()) { + LOG_INFO("Tracers: Group, 'Debug', exists PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: Group, 'Debug', does not exist FAIL"); + } + + I4 TotalLength = 0; + + for (std::string GroupName : GroupNames) { + std::pair GroupRange; + Err = Tracers::getGroupRange(GroupRange, GroupName); + + if (Err != 0) { + LOG_ERROR("Tracers: getGroupRange returns {} FAIL", Err); + RetVal += 1; + } + + auto [StartIndex, GroupLength] = GroupRange; + + TotalLength += GroupLength; + + // Check if a group contains more than one tracers + if (GroupLength > 0) { + LOG_INFO("Tracers: {} tracers retrieval PASS", GroupName); + } else { + RetVal += 1; + LOG_ERROR("Tracers: {} tracers retrieval FAIL", GroupName); + } + + // Check if tracer index is a member of the Group + for (I4 TracerIndex = StartIndex; + TracerIndex < StartIndex + GroupLength; ++TracerIndex) { + if (Tracers::isGroupMemberByIndex(TracerIndex, GroupName)) { + LOG_INFO("Tracers: {} group has the tracer index, {} PASS", + GroupName, TracerIndex); + } else { + RetVal += 1; + LOG_ERROR( + "Tracers: {} group does not have the tracer index, {} FAIL", + GroupName, TracerIndex); + } + } + + // Check if tracer index:name mapping is correct + for (I4 TracerIndex = StartIndex; + TracerIndex < StartIndex + GroupLength; ++TracerIndex) { + std::string TracerName; + Err = Tracers::getName(TracerName, TracerIndex); + if (Err != 0) { + LOG_ERROR("Tracers: getName returns {} FAIL", Err); + RetVal += 1; + } + + I4 RetTracerIndex; + Err = Tracers::getIndex(RetTracerIndex, TracerName); + if (Err != 0) { + LOG_ERROR("Tracers: getIndex returns {} FAIL", Err); + RetVal += 1; + } + + if (TracerIndex == RetTracerIndex) { + LOG_INFO("Tracers: {} group tracer:name mapping for {} is " + "correct PASS", + GroupName, TracerName); + } else { + RetVal += 1; + LOG_ERROR("Tracers: {} group tracer:name mapping for {} is not " + "correct FAIL", + GroupName, TracerName); + } + } + + // check if tracer has a corresponding field + for (I4 TracerIndex = StartIndex; + TracerIndex < StartIndex + GroupLength; ++TracerIndex) { + auto TracerField = Tracers::getFieldByIndex(TracerIndex); + + if (TracerField) { + LOG_INFO("Tracers: getFieldByIndex returns a field PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: getFieldByIndex returns nullptr FAIL"); + } + } + } + + I4 NTracers = Tracers::getNumTracers(); + I4 NCellsOwned = Tracers::NCellsOwned; + I4 NCellsAll = Tracers::NCellsAll; + I4 NCellsSize = Tracers::NCellsSize; + I4 NVertLevels = Tracers::NVertLevels; + + // Check if total number of tracers is correct + if (TotalLength == NTracers) { + LOG_INFO("Tracers: getNumTracers() returns correct tracer size PASS"); + } else { + RetVal += 1; + LOG_ERROR( + "Tracers: getNumTracers() returns incorrect tracer size FAIL"); + } + + // Reference host array of current time level + HostArray3DReal RefHostArray = + HostArray3DReal("RefHostArray", NTracers, NCellsSize, NVertLevels); + + // intialize tracer elements of all time levels + for (I4 TimeLevel = 0; TimeLevel + Tracers::NTimeLevels > 0; + --TimeLevel) { + HostArray3DReal TempHostArray; + Err = Tracers::getAllHost(TempHostArray, TimeLevel); + if (Err != 0) { + LOG_ERROR("getAllHost(TempHostArray, TimeLevel) returns non-zero " + "code: {}", + Err); + RetVal += 1; + } + + for (I4 Tracer = 0; Tracer < NTracers; ++Tracer) { + for (I4 Cell = 0; Cell < NCellsSize; Cell++) { + for (I4 Vert = 0; Vert < NVertLevels; Vert++) { + TempHostArray(Tracer, Cell, Vert) = + RefReal + Tracer + Cell + Vert + TimeLevel; + if (TimeLevel == 0) + RefHostArray(Tracer, Cell, Vert) = + TempHostArray(Tracer, Cell, Vert); + } + } + } + Tracers::copyToDevice(TimeLevel); + } + + // Reference device array of current time level + Array3DReal RefArray = + Array3DReal("RefArray", NTracers, NCellsSize, NVertLevels); + + Err = Tracers::getAll(RefArray, 0); + if (Err != 0) { + LOG_ERROR("getAll(RefArray, -1) returns non-zero code: {}", Err); + RetVal += 1; + } + + deepCopy(RefArray, RefArray); + + // Reference field data of all tracers + std::vector RefFieldDataArray; + + // get field references of all tracers + for (I4 Tracer = 0; Tracer < NTracers; ++Tracer) { + auto TracerField = Tracers::getFieldByIndex(Tracer); + RefFieldDataArray.push_back(TracerField->getDataArray()); + } + + // update time levels + Tracers::updateTimeLevels(); + + // getAll of previous time level(-1) should return the same to RefArray + Array3DReal PrevArray; + Err = Tracers::getAll(PrevArray, -1); + if (Err != 0) { + LOG_ERROR("getAll(PrevArray, -1) returns non-zero code: {}", Err); + RetVal += 1; + } + + count = -1; + + // check if time level shift works + parallelReduce( + "reduce1", {NTracers, NCellsOwned, NVertLevels}, + KOKKOS_LAMBDA(I4 Tracer, I4 Cell, I4 Vert, I4 & Accum) { + if (std::abs(PrevArray(Tracer, Cell, Vert) - + RefArray(Tracer, Cell, Vert)) > 1e-9) { + Accum++; + } + }, + count); + + if (count == 0) { + LOG_INFO("Tracers: Tracer data match after updateTimeLevels() PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: Not all tracer data match after " + "updateTimeLevels():{} FAIL", + count); + } + + // test getByName and getByIndex + for (I4 Tracer = 0; Tracer < NTracers; ++Tracer) { + std::string TracerName; + Tracers::getName(TracerName, Tracer); + + Array2DReal PrevTracer; + Err = Tracers::getByName(PrevTracer, -1, TracerName); + if (Err != 0) { + LOG_ERROR("getByName(PrevTracer, -1, TracerName) returns non-zero " + "code: {}", + Err); + RetVal += 1; + } + + count = -1; + + parallelReduce( + "reduce2", {NCellsOwned, NVertLevels}, + KOKKOS_LAMBDA(I4 Cell, I4 Vert, I4 & Accum) { + if (std::abs(PrevTracer(Cell, Vert) - + (RefReal + Tracer + Cell + Vert)) > 1e-9) { + Accum++; + } + }, + count); + + if (count == 0) { + LOG_INFO("Tracers: Tracer data from getByName match after " + "updateTimeLevels() PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: Not all tracer data from getByName match after " + "updateTimeLevels():{} FAIL", + count); + } + } + + // test field data - this test should generate positive count value. + // Referece and test field data are different due to updateTimeLevels() + for (I4 Tracer = 0; Tracer < NTracers; ++Tracer) { + + auto TracerField = Tracers::getFieldByIndex(Tracer); + Array2DReal TestFieldData = TracerField->getDataArray(); + Array2DReal RefFieldData = RefFieldDataArray[Tracer]; + + count = -1; + + parallelReduce( + "reduce3", {NCellsOwned, NVertLevels}, + KOKKOS_LAMBDA(I4 Cell, I4 Vert, I4 & Accum) { + if (std::abs(RefFieldData(Cell, Vert) - + TestFieldData(Cell, Vert)) > 1e-9) { + Accum++; + } + }, + count); + + if (count > 0) { + LOG_INFO("Tracers: Tracer field data correctly catch the " + "difference after updateTimeLevels() PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: Tracer field data should not match after " + "updateTimeLevels() FAIL"); + } + } + + // update time levels to cycle back to original index + for (I4 TimeLevel = -1; TimeLevel + Tracers::NTimeLevels > 0; + --TimeLevel) { + // update time levels + Tracers::updateTimeLevels(); + } + + // test field data - this test should generate zero count value + // Referece and test field data are the same because updateTimeLevels() is + // called NTimeLevels + for (I4 Tracer = 0; Tracer < NTracers; ++Tracer) { + auto TracerField = Tracers::getFieldByIndex(Tracer); + Array2DReal TestFieldData = TracerField->getDataArray(); + Array2DReal RefFieldData = RefFieldDataArray[Tracer]; + + count = -1; + + parallelReduce( + "reduce4", {NCellsOwned, NVertLevels}, + KOKKOS_LAMBDA(I4 Cell, I4 Vert, I4 & Accum) { + if (std::abs(RefFieldData(Cell, Vert) - + TestFieldData(Cell, Vert)) > 1e-9) { + Accum++; + } + }, + count); + + if (count == 0) { + LOG_INFO("Tracers: Tracer field data correctly match after " + "updateTimeLevels() back to original index PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: Not all tracer data match after " + "updateTimeLevels() back to original index FAIL"); + } + } + + count = 0; + + // Finally, check if getHostByName returns the correct host array + for (I4 Tracer = 0; Tracer < NTracers; ++Tracer) { + std::string TracerName; + Tracers::getName(TracerName, Tracer); + + HostArray2DReal TestHostArray; + Err = Tracers::getHostByName(TestHostArray, 0, TracerName); + if (Err != 0) { + LOG_ERROR("getHostByName(TestHostArray, 0, TracerName) returns " + "non-zero code: {}", + Err); + RetVal += 1; + } + + for (I4 Cell = 0; Cell < NCellsOwned; Cell++) { + for (I4 Vert = 0; Vert < NVertLevels; Vert++) { + if (std::abs(RefHostArray(Tracer, Cell, Vert) - + TestHostArray(Cell, Vert)) > 1e-9) + ++count; + } + } + } + + if (count == 0) { + LOG_INFO("Tracers: Tracer getHostByName correctly retreive tracer " + "data PASS"); + } else { + RetVal += 1; + LOG_ERROR("Tracers: Tracer getHostByName retreives incorrect tracer " + "data FAIL"); + } + + Tracers::clear(); + TimeStepper::clear(); + HorzMesh::clear(); + Decomp::clear(); + MachEnv::removeAll(); + FieldGroup::clear(); + Field::clear(); + Dimension::clear(); + + if (RetVal == 0) + LOG_INFO("Tracers: Successful completion"); + } + Kokkos::finalize(); + MPI_Finalize(); + + if (RetVal >= 256) + RetVal = 255; + + return RetVal; + +} // end of main +//===-----------------------------------------------------------------------===/