Skip to content

Commit

Permalink
Merge pull request #32 from MetOffice/feature/no_field_levels_calls
Browse files Browse the repository at this point in the history 
Feature/no field levels calls
  • Loading branch information
@phlndrwd
phlndrwd authored Jan 29, 2024
2 parents b7c95fa + a2f2a0a commit 65d6d7a
Show file tree
Hide file tree
Showing 5 changed files with 64 additions and 61 deletions.
25 changes: 13 additions & 12 deletions src/monio/AtlasReader.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,7 @@
#include "UtilsAtlas.h"
#include "Monio.h"

monio::AtlasReader::AtlasReader(const eckit::mpi::Comm& mpiCommunicator,
const int mpiRankOwner):
monio::AtlasReader::AtlasReader(const eckit::mpi::Comm& mpiCommunicator, const int mpiRankOwner):
mpiCommunicator_(mpiCommunicator),
mpiRankOwner_(mpiRankOwner) {
oops::Log::debug() << "AtlasReader::AtlasReader()" << std::endl;
Expand Down Expand Up @@ -116,14 +115,16 @@ void monio::AtlasReader::populateField(atlas::Field& field,
oops::Log::debug() << "AtlasReader::populateField()" << std::endl;
auto fieldView = atlas::array::make_view<T, 2>(field);
// Field with noFirstLevel == true should have been adjusted to have 70 levels.
if (noFirstLevel == true && field.levels() == consts::kVerticalFullSize) {
std::vector<atlas::idx_t> fieldShape = field.shape();
atlas::idx_t numLevels = fieldShape[consts::eVertical];
if (noFirstLevel == true && numLevels == consts::kVerticalFullSize) {
Monio::get().closeFiles();
utils::throwException("AtlasReader::populateField()> Field levels misconfiguration...");
// Only valid case for field with noFirstLevel == true. Field is adjusted to have 70 levels but
// read data still has enough to fill 71.
} else if (isLfricConvention == true &&
noFirstLevel == true &&
field.levels() == consts::kVerticalHalfSize) {
numLevels == consts::kVerticalHalfSize) {
for (int j = 1; j < consts::kVerticalFullSize; ++j) {
for (std::size_t i = 0; i < lfricToAtlasMap.size(); ++i) {
int index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size());
Expand All @@ -140,7 +141,7 @@ void monio::AtlasReader::populateField(atlas::Field& field,
// Valid case for isLfricConvention == false and fields noFirstLevel == false. Field is filled
// with all available data.
} else {
for (int j = 0; j < field.levels(); ++j) {
for (atlas::idx_t j = 0; j < numLevels; ++j) {
for (std::size_t i = 0; i < lfricToAtlasMap.size(); ++i) {
int index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size());
// Bounds checking
Expand Down Expand Up @@ -177,15 +178,14 @@ void monio::AtlasReader::populateField(atlas::Field& field,
const std::vector<T>& dataVec) {
oops::Log::debug() << "AtlasReader::populateField()" << std::endl;

std::vector<int> dimVec = field.shape();
std::vector<atlas::idx_t> fieldShape = field.shape();
if (field.metadata().get<bool>("global") == false) {
dimVec[consts::eHorizontal] = utilsatlas::getHorizontalSize(field);
fieldShape[consts::eHorizontal] = utilsatlas::getHorizontalSize(field);
}
auto fieldView = atlas::array::make_view<T, 2>(field);
int numLevels = field.levels();
for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) {
for (int j = 0; j < numLevels; ++j) {
int index = i + (j * dimVec[consts::eHorizontal]);
for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) {
for (atlas::idx_t j = 0; j < fieldShape[consts::eVertical]; ++j) {
atlas::idx_t index = i + (j * fieldShape[consts::eHorizontal]);
if (std::size_t(index) <= dataVec.size()) {
fieldView(i, j) = dataVec[index];
} else {
Expand All @@ -206,8 +206,9 @@ template void monio::AtlasReader::populateField<int>(atlas::Field& field,

atlas::Field monio::AtlasReader::getReadField(atlas::Field& field,
const bool noFirstLevel) {
std::vector<atlas::idx_t> fieldShape = field.shape();
// Check to ensure field has not been initialised with 71 levels
if (noFirstLevel == true && field.levels() == consts::kVerticalFullSize) {
if (noFirstLevel == true && fieldShape[consts::eVertical] == consts::kVerticalFullSize) {
atlas::array::DataType atlasType = field.datatype();
if (atlasType != atlasType.KIND_REAL64 &&
atlasType != atlasType.KIND_REAL32 &&
Expand Down
67 changes: 34 additions & 33 deletions src/monio/AtlasWriter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,7 @@
#include "UtilsAtlas.h"
#include "Writer.h"

monio::AtlasWriter::AtlasWriter(const eckit::mpi::Comm& mpiCommunicator,
const int mpiRankOwner):
monio::AtlasWriter::AtlasWriter(const eckit::mpi::Comm& mpiCommunicator, const int mpiRankOwner):
mpiCommunicator_(mpiCommunicator),
mpiRankOwner_(mpiRankOwner) {
oops::Log::debug() << "AtlasWriter::AtlasWriter()" << std::endl;
Expand Down Expand Up @@ -60,11 +59,11 @@ void monio::AtlasWriter::populateFileDataWithField(FileData& fileData,
Metadata& metadata = fileData.getMetadata();
Data& data = fileData.getData();
// Create dimensions
std::vector<int> dimVec = field.shape();
std::vector<atlas::idx_t> fieldShape = field.shape();
if (field.metadata().get<bool>("global") == false) {
dimVec[0] = utilsatlas::getHorizontalSize(field);
fieldShape[0] = utilsatlas::getHorizontalSize(field);
}
for (auto& dimSize : dimVec) {
for (auto& dimSize : fieldShape) {
std::string dimName = metadata.getDimensionName(dimSize);
if (dimName == consts::kNotFoundError) {
dimName = "dim" + std::to_string(dimCount_);
Expand All @@ -91,7 +90,7 @@ void monio::AtlasWriter::populateFileDataWithField(FileData& fileData,
metadata.addVariable(consts::kCoordVarNames[consts::eLongitude], lonVar);
metadata.addVariable(consts::kCoordVarNames[consts::eLatitude], latVar);

populateDataWithField(data, field, dimVec);
populateDataWithField(data, field, fieldShape);
addGlobalAttributes(metadata, false);
}
}
Expand Down Expand Up @@ -155,7 +154,7 @@ void monio::AtlasWriter::populateDataWithField(Data& data,

void monio::AtlasWriter::populateDataWithField(Data& data,
const atlas::Field& field,
const std::vector<int> dimensions) {
const std::vector<atlas::idx_t> dimensions) {
oops::Log::debug() << "AtlasWriter::populateDataWithField()" << std::endl;
std::shared_ptr<DataContainerBase> dataContainer = nullptr;
populateDataContainerWithField(dataContainer, field, dimensions);
Expand All @@ -170,7 +169,7 @@ void monio::AtlasWriter::populateDataContainerWithField(
oops::Log::debug() << "AtlasWriter::populateDataContainerWithField()" << std::endl;
if (mpiCommunicator_.rank() == mpiRankOwner_) {
atlas::array::DataType atlasType = field.datatype();
int fieldSize = utilsatlas::getGlobalDataSize(field);
atlas::idx_t fieldSize = utilsatlas::getGlobalDataSize(field);
switch (atlasType.kind()) {
case atlasType.KIND_INT32: {
if (dataContainer == nullptr) {
Expand Down Expand Up @@ -217,13 +216,13 @@ void monio::AtlasWriter::populateDataContainerWithField(
void monio::AtlasWriter::populateDataContainerWithField(
std::shared_ptr<monio::DataContainerBase>& dataContainer,
const atlas::Field& field,
const std::vector<int>& dimensions) {
const std::vector<atlas::idx_t>& dimensions) {
oops::Log::debug() << "AtlasWriter::populateDataContainerWithField()" << std::endl;
oops::Log::info() << "AtlasWriter::populateDataContainerWithField()" << std::endl;
if (mpiCommunicator_.rank() == mpiRankOwner_) {
std::string fieldName = field.name();
atlas::array::DataType atlasType = field.datatype();
int fieldSize = utilsatlas::getGlobalDataSize(field);
atlas::idx_t fieldSize = utilsatlas::getGlobalDataSize(field);
switch (atlasType.kind()) {
case atlasType.KIND_INT32: {
if (dataContainer == nullptr) {
Expand Down Expand Up @@ -272,16 +271,17 @@ void monio::AtlasWriter::populateDataVec(std::vector<T>& dataVec,
const atlas::Field& field,
const std::vector<size_t>& lfricToAtlasMap) {
oops::Log::debug() << "AtlasWriter::populateDataVec() " << field.name() << std::endl;
int numLevels = field.levels();
std::vector<atlas::idx_t> fieldShape = field.shape();
atlas::idx_t numLevels = fieldShape[consts::eVertical];
if ((lfricToAtlasMap.size() * numLevels) != dataVec.size()) {
Monio::get().closeFiles();
utils::throwException("AtlasWriter::populateDataVec()> "
"Data container is not configured for the expected data...");
}
auto fieldView = atlas::array::make_view<T, 2>(field);
for (std::size_t i = 0; i < lfricToAtlasMap.size(); ++i) {
for (int j = 0; j < numLevels; ++j) {
int index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size());
for (atlas::idx_t j = 0; j < numLevels; ++j) {
atlas::idx_t index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size());
dataVec[index] = fieldView(i, j);
}
}
Expand All @@ -300,26 +300,26 @@ template void monio::AtlasWriter::populateDataVec<int>(std::vector<int>& dataVec
template<typename T>
void monio::AtlasWriter::populateDataVec(std::vector<T>& dataVec,
const atlas::Field& field,
const std::vector<int>& dimensions) {
const std::vector<atlas::idx_t>& dimensions) {
oops::Log::debug() << "AtlasWriter::populateDataVec()" << std::endl;
auto fieldView = atlas::array::make_view<T, 2>(field);
for (int i = 0; i < dimensions[consts::eHorizontal]; ++i) {
for (int j = 0; j < dimensions[consts::eVertical]; ++j) {
int index = j + (i * dimensions[consts::eVertical]);
for (atlas::idx_t i = 0; i < dimensions[consts::eHorizontal]; ++i) {
for (atlas::idx_t j = 0; j < dimensions[consts::eVertical]; ++j) {
atlas::idx_t index = j + (i * dimensions[consts::eVertical]);
dataVec[index] = fieldView(i, j);
}
}
}

template void monio::AtlasWriter::populateDataVec<double>(std::vector<double>& dataVec,
const atlas::Field& field,
const std::vector<int>& dimensions);
const std::vector<atlas::idx_t>& dimensions);
template void monio::AtlasWriter::populateDataVec<float>(std::vector<float>& dataVec,
const atlas::Field& field,
const std::vector<int>& dimensions);
const std::vector<atlas::idx_t>& dimensions);
template void monio::AtlasWriter::populateDataVec<int>(std::vector<int>& dataVec,
const atlas::Field& field,
const std::vector<int>& dimensions);
const std::vector<atlas::idx_t>& dimensions);

atlas::Field monio::AtlasWriter::getWriteField(atlas::Field& field,
const std::string& writeName,
Expand All @@ -333,14 +333,15 @@ atlas::Field monio::AtlasWriter::getWriteField(atlas::Field& field,
Monio::get().closeFiles();
utils::throwException("AtlasWriter::getWriteField())> Data type not coded for...");
}
std::vector<atlas::idx_t> fieldShape = field.shape();
// Erroneous case. For noFirstLevel == true field should have 70 levels
if (noFirstLevel == true && field.levels() == consts::kVerticalFullSize) {
if (noFirstLevel == true && fieldShape[consts::eVertical] == consts::kVerticalFullSize) {
Monio::get().closeFiles();
utils::throwException("AtlasWriter::getWriteField()> Field levels misconfiguration...");
}
// WARNING - This name-check is an LFRic-Lite specific convention...
if (utils::findInVector(consts::kMissingVariableNames, writeName) == false) {
if (noFirstLevel == true && field.levels() == consts::kVerticalHalfSize) {
if (noFirstLevel == true && fieldShape[consts::eVertical] == consts::kVerticalHalfSize) {
atlas::util::Config atlasOptions = atlas::option::name(writeName) |
atlas::option::global(0) |
atlas::option::levels(consts::kVerticalFullSize);
Expand Down Expand Up @@ -373,14 +374,14 @@ atlas::Field monio::AtlasWriter::copySurfaceLevel(const atlas::Field& inputField
atlas::Field copiedField = functionSpace.createField<T>(atlasOptions);
auto copiedFieldView = atlas::array::make_view<T, 2>(copiedField);
auto inputFieldView = atlas::array::make_view<T, 2>(inputField);
std::vector<int> dimVec = inputField.shape();
for (int j = 0; j < dimVec[consts::eVertical]; ++j) {
for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) {
std::vector<atlas::idx_t> fieldShape = inputField.shape();
for (atlas::idx_t j = 0; j < fieldShape[consts::eVertical]; ++j) {
for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) {
copiedFieldView(i, j + 1) = inputFieldView(i, j);
}
}
// Copy surface level of input field
for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) {
for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) {
copiedFieldView(i, 0) = inputFieldView(i, 0);
}
return copiedField;
Expand All @@ -400,16 +401,16 @@ void monio::AtlasWriter::addVariableDimensions(const atlas::Field& field,
const Metadata& metadata,
std::shared_ptr<monio::Variable> var,
const std::string& vertConfigName) {
std::vector<int> dimVec = field.shape();
if (field.metadata().get<bool>("global") == false) {
dimVec[0] = utilsatlas::getHorizontalSize(field); // If so, get the 2D size of the Field
std::vector<atlas::idx_t> fieldShape = field.shape();
if (field.metadata().get<bool>("global") == false) { // If so, get the 2D size of the Field
fieldShape[consts::eHorizontal] = utilsatlas::getHorizontalSize(field);
}
// Reversal of dims required for LFRic files. Currently applied to all output files.
std::reverse(dimVec.begin(), dimVec.end());
for (auto& dimSize : dimVec) {
std::reverse(fieldShape.begin(), fieldShape.end());
for (auto& dimSize : fieldShape) {
std::string dimName;
if (vertConfigName != "" && metadata.isDimDefined(vertConfigName) &&
dimSize == metadata.getDimension(vertConfigName)) {
dimSize == metadata.getDimension(vertConfigName)) {
dimName = vertConfigName;
} else {
dimName = metadata.getDimensionName(dimSize);
Expand All @@ -424,7 +425,7 @@ void monio::AtlasWriter::addGlobalAttributes(Metadata& metadata, const bool isLf
// Initialise variables
std::string variableConvention =
isLfricConvention == true ? consts::kNamingConventions[consts::eLfricConvention] :
consts::kNamingConventions[consts::eJediConvention];
consts::kNamingConventions[consts::eJediConvention];
// Create attribute objects
std::shared_ptr<monio::AttributeString> namingAttr =
std::make_shared<AttributeString>(std::string(consts::kVariableConventionName),
Expand Down
4 changes: 2 additions & 2 deletions src/monio/AtlasWriter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ class AtlasWriter {
/// populateFileDataWithField where metadata are created.
void populateDataWithField(Data& data,
const atlas::Field& field,
const std::vector<int> dimensions);
const std::vector<atlas::idx_t> dimensions);

/// \brief Derives the container type and makes the call to populate it. Used where metadata are
/// provided and data are written in LFRic order.
Expand All @@ -101,7 +101,7 @@ class AtlasWriter {
/// \brief Iterates through field and populates vector with data from field in Atlas order.
template<typename T> void populateDataVec(std::vector<T>& dataVec,
const atlas::Field& field,
const std::vector<int>& dimensions);
const std::vector<atlas::idx_t>& dimensions);

/// \brief Map JEDI fields back into LFRic function space.
atlas::Field getWriteField(atlas::Field& inputField,
Expand Down
25 changes: 13 additions & 12 deletions src/monio/UtilsAtlas.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ std::vector<atlas::PointLonLat> getAtlasCoords(const atlas::Field& field) {
if (field.metadata().get<bool>("global") == false) {
auto lonLatField = field.functionspace().lonlat();
auto lonLatView = atlas::array::make_view<double, 2>(lonLatField);
for (int i = 0; i < getHorizontalSize(field); ++i) {
for (atlas::idx_t i = 0; i < getHorizontalSize(field); ++i) {
atlasCoords.push_back(atlas::PointLonLat(lonLatView(i, consts::eLongitude),
lonLatView(i, consts::eLatitude)));
}
Expand Down Expand Up @@ -131,8 +131,9 @@ std::vector<size_t> createLfricAtlasMap(const std::vector<atlas::PointLonLat>& a
atlas::Field getGlobalField(const atlas::Field& field) {
if (field.metadata().get<bool>("global") == false) {
atlas::array::DataType atlasType = field.datatype();
std::vector<atlas::idx_t> fieldShape = field.shape();
atlas::util::Config atlasOptions = atlas::option::name(field.name()) |
atlas::option::levels(field.levels()) |
atlas::option::levels(fieldShape[consts::eVertical]) |
atlas::option::datatype(atlasType) |
atlas::option::global(0);
if (atlasType != atlasType.KIND_REAL64 &&
Expand Down Expand Up @@ -164,11 +165,11 @@ atlas::FieldSet getGlobalFieldSet(const atlas::FieldSet& fieldSet) {
}
}

int getHorizontalSize(const atlas::Field& field) {
atlas::idx_t getHorizontalSize(const atlas::Field& field) {
atlas::Field ghostField = field.functionspace().ghost();
int size = 0;
atlas::idx_t size = 0;
auto ghostView = atlas::array::make_view<int, 1>(ghostField);
for (int i = ghostField.size() - 1; i >= 0; --i) {
for (atlas::idx_t i = ghostField.size() - 1; i >= 0; --i) {
if (ghostView(i) == 0) {
size = i + 1;
break;
Expand All @@ -177,10 +178,10 @@ int getHorizontalSize(const atlas::Field& field) {
return size;
}

int getGlobalDataSize(const atlas::Field& field) {
std::vector<int> dimVec = field.shape();
int size = 1;
for (const auto& dim : dimVec) {
atlas::idx_t getGlobalDataSize(const atlas::Field& field) {
std::vector<atlas::idx_t> fieldShape = field.shape();
atlas::idx_t size = 1;
for (const auto& dim : fieldShape) {
size *= dim;
}
return size;
Expand Down Expand Up @@ -216,9 +217,9 @@ bool compareFieldSets(const atlas::FieldSet& aSet, const atlas::FieldSet& bSet)
bool compareFields(const atlas::Field& a, const atlas::Field& b) {
const auto aView = atlas::array::make_view<const double, 2>(a);
const auto bView = atlas::array::make_view<const double, 2>(b);
std::vector<int> dimVec = a.shape();
for (int j = 0; j < dimVec[consts::eVertical]; ++j) {
for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) {
std::vector<atlas::idx_t> fieldShape = a.shape();
for (atlas::idx_t j = 0; j < fieldShape[consts::eVertical]; ++j) {
for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) {
if (aView(i, j) != bView(i, j)) {
return false;
}
Expand Down
4 changes: 2 additions & 2 deletions src/monio/UtilsAtlas.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,8 +43,8 @@ namespace utilsatlas {

atlas::Field getGlobalField(const atlas::Field& field);

int getHorizontalSize(const atlas::Field& field); // Just 2D size. Any field.
int getGlobalDataSize(const atlas::Field& field); // Full 3D size of data. Global fields only
atlas::idx_t getHorizontalSize(const atlas::Field& field); // Just 2D size. Any field.
atlas::idx_t getGlobalDataSize(const atlas::Field& field); // Full 3D size of global field.

int atlasTypeToMonioEnum(atlas::array::DataType atlasType);

Expand Down

0 comments on commit 65d6d7a

Please sign in to comment.