Skip to content

Enable choice of LFRic-Atlas map #54

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
Chris Thomas (ctgh) wants to merge 1 commit into from
Closed

Enable choice of LFRic-Atlas map #54

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions src/monio/Constants.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,11 @@ enum eDataTypes {
eNumberOfDataTypes
};

enum eLfricAtlasMapMethods {
eTrivial,
eKDtree
};

/// String/Views ///////////////////////////////////////////////////////////////////////////////////

const std::string_view kTimeDimName = "time_counter";
Expand Down
28 changes: 18 additions & 10 deletions src/monio/Monio.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,10 @@ monio::Monio::~Monio() {
}

void monio::Monio::readState(atlas::FieldSet& localFieldSet,
const std::vector<consts::FieldMetadata>& fieldMetadataVec,
const std::string& filePath,
const util::DateTime& dateTime) {
const std::vector<consts::FieldMetadata>& fieldMetadataVec,
const std::string& filePath,
const util::DateTime& dateTime,
const consts::eLfricAtlasMapMethods& lfricAtlasMapMethod) {
oops::Log::trace() << "Monio::readState()" << std::endl;
if (localFieldSet.size() == 0) {
Monio::get().closeFiles();
Expand All @@ -62,7 +63,8 @@ void monio::Monio::readState(atlas::FieldSet& localFieldSet,
auto& functionSpace = globalField.functionspace();
auto grid = utilsatlas::getGridFromFunctionSpace(functionSpace);
// Initialise file
int variableConvention = initialiseFile(grid, filePath, true);
const int variableConvention = initialiseFile(grid, filePath, true,
lfricAtlasMapMethod);
// getFileData returns a copy of FileData (with required LFRic mesh data), so read data
// is discarded when FileData goes out-of-scope for reading subsequent fields.
FileData fileData = getFileData(grid.name());
Expand Down Expand Up @@ -106,7 +108,8 @@ void monio::Monio::readState(atlas::FieldSet& localFieldSet,

void monio::Monio::readIncrements(atlas::FieldSet& localFieldSet,
const std::vector<consts::FieldMetadata>& fieldMetadataVec,
const std::string& filePath) {
const std::string& filePath,
const consts::eLfricAtlasMapMethods& lfricAtlasMapMethod) {
oops::Log::trace() << "Monio::readIncrements()" << std::endl;
if (localFieldSet.size() == 0) {
Monio::get().closeFiles();
Expand All @@ -123,7 +126,8 @@ void monio::Monio::readIncrements(atlas::FieldSet& localFieldSet,
auto grid = utilsatlas::getGridFromFunctionSpace(functionSpace);

// Initialise file
int variableConvention = initialiseFile(grid, filePath);
const int variableConvention = initialiseFile(grid, filePath, false,
lfricAtlasMapMethod);
// getFileData returns a copy of FileData (with required LFRic mesh data), so read data
// is discarded when FileData goes out-of-scope for reading subsequent fields.
FileData fileData = getFileData(grid.name());
Expand Down Expand Up @@ -329,7 +333,8 @@ void monio::Monio::closeFiles() {

int monio::Monio::initialiseFile(const atlas::Grid& grid,
const std::string& filePath,
bool doCreateDateTimes) {
const bool doCreateDateTimes,
const consts::eLfricAtlasMapMethods& lfricAtlasMapMethod) {
oops::Log::trace() << "Monio::initialiseFile()" << std::endl;
int variableConvention = consts::eLfricConvention; // LFRic convention is default
if (mpiCommunicator_.rank() == mpiRankOwner_) {
Expand All @@ -343,7 +348,7 @@ int monio::Monio::initialiseFile(const atlas::Grid& grid,
reader_.readFullDatum(fileData, std::string(consts::kVerticalFullName));
reader_.readFullDatum(fileData, std::string(consts::kVerticalHalfName));
// Process read data
createLfricAtlasMap(fileData, grid);
createLfricAtlasMap(fileData, grid, lfricAtlasMapMethod);
if (doCreateDateTimes == true) {
reader_.readFullDatum(fileData, std::string(consts::kTimeVarName));
createDateTimes(fileData,
Expand Down Expand Up @@ -390,7 +395,9 @@ monio::FileData monio::Monio::getFileData(const std::string& gridName) {
return FileData(); // This function is called by all PEs. A return is essential.
}

void monio::Monio::createLfricAtlasMap(FileData& fileData, const atlas::Grid& grid) {
void monio::Monio::createLfricAtlasMap(FileData& fileData,
const atlas::Grid& grid,
const consts::eLfricAtlasMapMethods& lfricAtlasMapMethod) {
oops::Log::trace() << "Monio::createLfricAtlasMap()" << std::endl;
if (mpiCommunicator_.rank() == mpiRankOwner_) {
if (fileData.getLfricAtlasMap().size() == 0) {
Expand All @@ -399,7 +406,8 @@ void monio::Monio::createLfricAtlasMap(FileData& fileData, const atlas::Grid& gr
reader_.getCoordData(fileData, consts::kLfricCoordVarNames);
std::vector<atlas::PointLonLat> lfricCoords = utilsatlas::getLfricCoords(coordData);
std::vector<atlas::PointLonLat> atlasCoords = utilsatlas::getAtlasCoords(grid);
fileData.setLfricAtlasMap(utilsatlas::createLfricAtlasMap(atlasCoords, lfricCoords));
fileData.setLfricAtlasMap(utilsatlas::createLfricAtlasMap(atlasCoords, lfricCoords,
lfricAtlasMapMethod));
}
}
}
Expand Down
18 changes: 14 additions & 4 deletions src/monio/Monio.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@

#include "AtlasReader.h"
#include "AtlasWriter.h"
#include "Constants.h"
#include "FileData.h"
#include "Reader.h"
#include "Writer.h"
Expand All @@ -41,12 +42,16 @@ class Monio {
void readState(atlas::FieldSet& localFieldSet,
const std::vector<consts::FieldMetadata>& fieldMetadataVec,
const std::string& filePath,
const util::DateTime& dateTime);
const util::DateTime& dateTime,
const consts::eLfricAtlasMapMethods& lfricAtlasMap =
consts::eLfricAtlasMapMethods::eKDtree);

/// \brief Reads files without a time component, i.e. increment files.
void readIncrements(atlas::FieldSet& localFieldSet,
const std::vector<consts::FieldMetadata>& fieldMetadataVec,
const std::string& filePath);
const std::string& filePath,
const consts::eLfricAtlasMapMethods& lfricAtlasMap =
consts::eLfricAtlasMapMethods::eKDtree);

/// \brief Writes increment files. No time component but the variables can use JEDI or LFRic write
/// names.
Expand All @@ -73,7 +78,9 @@ class Monio {
/// it's called from LFRic-Lite.
int initialiseFile(const atlas::Grid& grid,
const std::string& filePath,
bool doCreateDateTimes = false);
const bool doCreateDateTimes = false,
const consts::eLfricAtlasMapMethods& lfricAtlasMap =
consts::eLfricAtlasMapMethods::eKDtree);

private:
/// \brief Private class constructor to prevent instantiation outside of the singleton.
Expand All @@ -89,7 +96,10 @@ class Monio {
FileData getFileData(const std::string& gridName);

/// \brief Creates and stores a map between Atlas and LFRic horizontal ordering.
void createLfricAtlasMap(FileData& fileData, const atlas::Grid& grid);
void createLfricAtlasMap(FileData& fileData,
const atlas::Grid& grid,
const consts::eLfricAtlasMapMethods& lfricAtlasMap =
consts::eLfricAtlasMapMethods::eKDtree);

/// \brief Creates and stores date-times from a state file.
void createDateTimes(FileData& fileData,
Expand Down
33 changes: 19 additions & 14 deletions src/monio/UtilsAtlas.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,8 @@ std::vector<std::shared_ptr<monio::DataContainerBase>> convertLatLonToContainers
}

std::vector<size_t> createLfricAtlasMap(const std::vector<atlas::PointLonLat>& atlasCoords,
const std::vector<atlas::PointLonLat>& lfricCoords) {
const std::vector<atlas::PointLonLat>& lfricCoords,
const consts::eLfricAtlasMapMethods& lfricAtlasMapMethod) {
std::vector<size_t> lfricAtlasMap;
// Essential check to ensure grid is configured to accommodate the data
if (atlasCoords.size() != lfricCoords.size()) {
Expand All @@ -108,22 +109,26 @@ std::vector<size_t> createLfricAtlasMap(const std::vector<atlas::PointLonLat>& a
}
lfricAtlasMap.reserve(atlasCoords.size());

// Make a kd-tree using atlasLonLat as the point,
// with element index i as payload
std::vector<size_t> indices(atlasCoords.size());
std::iota(begin(indices), end(indices), 0);

const atlas::Geometry unitSphere(1.0);
atlas::util::IndexKDTree tree(unitSphere);
tree.build(atlasCoords, indices);

// Partitioning vector to create atlas::Distribution for mesh generation
std::vector<int> partitioning(tree.size(), -1);

// find atlas global indices for each element of modelLonLat
for (const auto& lfricCoord : lfricCoords) {
auto idx = tree.closestPoint(lfricCoord).payload();
lfricAtlasMap.push_back(idx);
// Trivially map between the LFRic and Atlas indices.
if (lfricAtlasMapMethod == consts::eLfricAtlasMapMethods::eTrivial) {
lfricAtlasMap = indices;
} else if (lfricAtlasMapMethod == consts::eLfricAtlasMapMethods::eKDtree) {
Comment on lines +116 to +118
Copy link
Copy Markdown
Contributor

@mo-joshuacolclough Josh Colclough (mo-joshuacolclough) Jul 7, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice to see a switch (lfricAtlasMapMethod) { case ... } for this (very minor change!) - but maybe worth holding off until we understand the issue

// Make a kd-tree using atlasLonLat as the point,
// with element index i as payload
const atlas::Geometry unitSphere(1.0);
atlas::util::IndexKDTree tree(unitSphere);
tree.build(atlasCoords, indices);

// find atlas global indices for each element of modelLonLat
for (const auto& lfricCoord : lfricCoords) {
const auto idx = tree.closestPoint(lfricCoord).payload();
lfricAtlasMap.push_back(idx);
}
} else {
utils::throwException("utilsatlas::createLfricAtlasMap()> Invalid LFRic-Atlas map chosen...");
}
return lfricAtlasMap;
}
Expand Down
3 changes: 2 additions & 1 deletion src/monio/UtilsAtlas.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,8 @@ namespace utilsatlas {
const std::vector<std::string>& coordNames);

std::vector<size_t> createLfricAtlasMap(const std::vector<atlas::PointLonLat>& atlasCoords,
const std::vector<atlas::PointLonLat>& lfricCoords);
const std::vector<atlas::PointLonLat>& lfricCoords,
const consts::eLfricAtlasMapMethods& lfricAtlasMapMethod);

atlas::FieldSet getGlobalFieldSet(const atlas::FieldSet& fieldSet);

Expand Down