|
moab
|
Child helper class for Eulerian Spectral grid (CAM_EUL) More...
#include <NCHelperEuler.hpp>
Public Member Functions | |
| NCHelperEuler (ReadNC *readNC, int fileId, const FileOptions &opts, EntityHandle fileSet) | |
Static Public Member Functions | |
| static bool | can_read_file (ReadNC *readNC, int fileId) |
Private Member Functions | |
| virtual ErrorCode | init_mesh_vals () |
| Interfaces to be implemented in child classes. | |
| virtual std::string | get_mesh_type_name () |
Child helper class for Eulerian Spectral grid (CAM_EUL)
Definition at line 17 of file NCHelperEuler.hpp.
| moab::NCHelperEuler::NCHelperEuler | ( | ReadNC * | readNC, |
| int | fileId, | ||
| const FileOptions & | opts, | ||
| EntityHandle | fileSet | ||
| ) | [inline] |
Definition at line 20 of file NCHelperEuler.hpp.
: ScdNCHelper(readNC, fileId, opts, fileSet) {}
| bool moab::NCHelperEuler::can_read_file | ( | ReadNC * | readNC, |
| int | fileId | ||
| ) | [static] |
Definition at line 16 of file NCHelperEuler.cpp.
{
std::vector<std::string>& dimNames = readNC->dimNames;
// If dimension names "lon" AND "lat' exist then it could be either the Eulerian Spectral grid or the FV grid
if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
dimNames.end(), std::string("lat")) != dimNames.end())) {
// If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV grid
if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
dimNames.end(), std::string("slat")) != dimNames.end()))
return false;
// Make sure it is CAM grid
std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
if (attIt == readNC->globalAtts.end()) {
readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
return false;
}
unsigned int sz = attIt->second.attLen;
std::string att_data;
att_data.resize(sz + 1);
att_data[sz] = '\000';
int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
if (success != 0) {
readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
return false;
}
if (att_data.find("CAM") == std::string::npos)
return false;
return true;
}
return false;
}
| virtual std::string moab::NCHelperEuler::get_mesh_type_name | ( | ) | [inline, private, virtual] |
| ErrorCode moab::NCHelperEuler::init_mesh_vals | ( | ) | [private, virtual] |
Interfaces to be implemented in child classes.
Implements moab::NCHelper.
Definition at line 52 of file NCHelperEuler.cpp.
{
Interface*& mbImpl = _readNC->mbImpl;
std::vector<std::string>& dimNames = _readNC->dimNames;
std::vector<int>& dimLens = _readNC->dimLens;
std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
DebugOutput& dbgOut = _readNC->dbgOut;
bool& isParallel = _readNC->isParallel;
int& partMethod = _readNC->partMethod;
ScdParData& parData = _readNC->parData;
// Look for names of center i/j dimensions
// First i
std::vector<std::string>::iterator vit;
unsigned int idx;
if ((vit = std::find(dimNames.begin(), dimNames.end(), "lon")) != dimNames.end())
idx = vit - dimNames.begin();
else {
ERRORR(MB_FAILURE, "Couldn't find 'lon' dimension.");
}
iCDim = idx;
gCDims[0] = 0;
gCDims[3] = dimLens[idx] - 1;
// Check i periodicity and set globallyPeriodic[0]
std::vector<double> til_vals(2);
ErrorCode rval = read_coordinate("lon", gCDims[3] - 1, gCDims[3], til_vals);
ERRORR(rval, "Trouble reading 'lon' variable.");
if (std::fabs(2 * til_vals[1] - til_vals[0] - 360) < 0.001)
globallyPeriodic[0] = 1;
// Now we can set gDims for i
gDims[0] = 0;
gDims[3] = gCDims[3] + (globallyPeriodic[0] ? 0 : 1); // Only if not periodic is vertex param max > elem param max
// Then j
if ((vit = std::find(dimNames.begin(), dimNames.end(), "lat")) != dimNames.end())
idx = vit - dimNames.begin();
else {
ERRORR(MB_FAILURE, "Couldn't find 'lat' dimension.");
}
jCDim = idx;
gCDims[1] = 0;
gCDims[4] = dimLens[idx] - 1;
// For Eul models, will always be non-periodic in j
gDims[1] = 0;
gDims[4] = gCDims[4] + 1;
// Try a truly 2D mesh
gDims[2] = -1;
gDims[5] = -1;
// Look for time dimension
if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
idx = vit - dimNames.begin();
else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
idx = vit - dimNames.begin();
else {
ERRORR(MB_FAILURE, "Couldn't find 'time' or 't' dimension.");
}
tDim = idx;
nTimeSteps = dimLens[idx];
// Look for level dimension
if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
idx = vit - dimNames.begin();
else if ((vit = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end())
idx = vit - dimNames.begin();
else {
ERRORR(MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension.");
}
levDim = idx;
nLevels = dimLens[idx];
// Parse options to get subset
int rank = 0, procs = 1;
#ifdef USE_MPI
if (isParallel) {
ParallelComm*& myPcomm = _readNC->myPcomm;
rank = myPcomm->proc_config().proc_rank();
procs = myPcomm->proc_config().proc_size();
}
#endif
if (procs > 1) {
for (int i = 0; i < 6; i++)
parData.gDims[i] = gDims[i];
for (int i = 0; i < 3; i++)
parData.gPeriodic[i] = globallyPeriodic[i];
parData.partMethod = partMethod;
int pdims[3];
rval = ScdInterface::compute_partition(procs, rank, parData, lDims, locallyPeriodic, pdims);
if (MB_SUCCESS != rval)
return rval;
for (int i = 0; i < 3; i++)
parData.pDims[i] = pdims[i];
dbgOut.tprintf(1, "Partition: %dx%d (out of %dx%d)\n",
lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1);
if (0 == rank)
dbgOut.tprintf(1, "Contiguous chunks of size %d bytes.\n", 8 * (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1));
}
else {
for (int i = 0; i < 6; i++)
lDims[i] = gDims[i];
locallyPeriodic[0] = globallyPeriodic[0];
}
_opts.get_int_option("IMIN", lDims[0]);
_opts.get_int_option("IMAX", lDims[3]);
_opts.get_int_option("JMIN", lDims[1]);
_opts.get_int_option("JMAX", lDims[4]);
// Now get actual coordinate values for vertices and cell centers
lCDims[0] = lDims[0];
if (locallyPeriodic[0])
// If locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem coords
lCDims[3] = lDims[3];
else
lCDims[3] = lDims[3] - 1;
// For Eul models, will always be non-periodic in j
lCDims[1] = lDims[1];
lCDims[4] = lDims[4] - 1;
// Resize vectors to store values later
if (-1 != lDims[0])
ilVals.resize(lDims[3] - lDims[0] + 1);
if (-1 != lCDims[0])
ilCVals.resize(lCDims[3] - lCDims[0] + 1);
if (-1 != lDims[1])
jlVals.resize(lDims[4] - lDims[1] + 1);
if (-1 != lCDims[1])
jlCVals.resize(lCDims[4] - lCDims[1] + 1);
if (nTimeSteps > 0)
tVals.resize(nTimeSteps);
// Now read coord values
std::map<std::string, ReadNC::VarData>::iterator vmit;
if (-1 != lCDims[0]) {
if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
rval = read_coordinate("lon", lCDims[0], lCDims[3], ilCVals);
ERRORR(rval, "Trouble reading 'lon' variable.");
}
else {
ERRORR(MB_FAILURE, "Couldn't find 'lon' variable.");
}
}
if (-1 != lCDims[1]) {
if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
rval = read_coordinate("lat", lCDims[1], lCDims[4], jlCVals);
ERRORR(rval, "Trouble reading 'lat' variable.");
}
else {
ERRORR(MB_FAILURE, "Couldn't find 'lat' variable.");
}
}
if (-1 != lDims[0]) {
if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
double dif = (ilCVals[1] - ilCVals[0]) / 2;
std::size_t i;
for (i = 0; i != ilCVals.size(); i++)
ilVals[i] = ilCVals[i] - dif;
// The last one is needed only if not periodic
if (!locallyPeriodic[0])
ilVals[i] = ilCVals[i - 1] + dif;
}
else {
ERRORR(MB_FAILURE, "Couldn't find 'lon' variable.");
}
}
if (-1 != lDims[1]) {
if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
if (!isParallel || ((gDims[4] - gDims[1]) == (lDims[4] - lDims[1]))) {
std::string gwName("gw");
std::vector<double> gwVals(lDims[4] - lDims[1] - 1);
rval = read_coordinate(gwName.c_str(), lDims[1], lDims[4] - 2, gwVals);
ERRORR(rval, "Trouble reading 'gw' variable.");
// Copy the correct piece
jlVals[0] = -(M_PI / 2) * 180 / M_PI;
std::size_t i = 0;
double gwSum = -1;
for (i = 1; i != gwVals.size() + 1; i++) {
gwSum += gwVals[i - 1];
jlVals[i] = std::asin(gwSum) * 180 / M_PI;
}
jlVals[i] = 90.0; // Using value of i after loop exits.
}
else {
std::string gwName("gw");
double gwSum = 0;
// If this is the first row
if (lDims[1] == gDims[1]) {
std::vector<double> gwVals(lDims[4]);
rval = read_coordinate(gwName.c_str(), 0, lDims[4] - 1, gwVals);
ERRORR(rval, "Trouble reading 'gw' variable.");
// Copy the correct piece
jlVals[0] = -(M_PI / 2) * 180 / M_PI;
gwSum = -1;
for (std::size_t i = 1; i != jlVals.size(); i++) {
gwSum += gwVals[i - 1];
jlVals[i] = std::asin(gwSum) * 180 / M_PI;
}
}
// Or if it's the last row
else if (lDims[4] == gDims[4]) {
std::vector<double> gwVals(lDims[4] - 1);
rval = read_coordinate(gwName.c_str(), 0, lDims[4] - 2, gwVals);
ERRORR(rval, "Trouble reading 'gw' variable.");
// copy the correct piece
gwSum = -1;
for (int j = 0; j != lDims[1] - 1; j++)
gwSum += gwVals[j];
std::size_t i = 0;
for (; i != jlVals.size() - 1; i++) {
gwSum += gwVals[lDims[1] - 1 + i];
jlVals[i] = std::asin(gwSum) * 180 / M_PI;
}
jlVals[i] = 90.0; // Using value of i after loop exits.
}
// It's in the middle
else {
int start = lDims[1] - 1;
int end = lDims[4] - 1;
std::vector<double> gwVals(end);
rval = read_coordinate(gwName.c_str(), 0, end - 1, gwVals);
ERRORR(rval, "Trouble reading 'gw' variable.");
gwSum = -1;
for (int j = 0; j != start - 1; j++)
gwSum += gwVals[j];
std::size_t i = 0;
for (; i != jlVals.size(); i++) {
gwSum += gwVals[start - 1 + i];
jlVals[i] = std::asin(gwSum) * 180 / M_PI;
}
}
}
}
else {
ERRORR(MB_FAILURE, "Couldn't find 'lat' variable.");
}
}
// Store time coordinate values in tVals
if (nTimeSteps > 0) {
if ((vmit = varInfo.find("time")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
rval = read_coordinate("time", 0, nTimeSteps - 1, tVals);
ERRORR(rval, "Trouble reading 'time' variable.");
}
else if ((vmit = varInfo.find("t")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
rval = read_coordinate("t", 0, nTimeSteps - 1, tVals);
ERRORR(rval, "Trouble reading 't' variable.");
}
else {
// If expected time variable is not available, set dummy time coordinate values to tVals
for (int t = 0; t < nTimeSteps; t++)
tVals.push_back((double)t);
}
}
dbgOut.tprintf(1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4]);
dbgOut.tprintf(1, "%d elements, %d vertices\n", (lDims[3] - lDims[0]) * (lDims[4] - lDims[1]), (lDims[3] - lDims[0] + 1)
* (lDims[4] - lDims[1] + 1));
// For each variable, determine the entity location type and number of levels
std::map<std::string, ReadNC::VarData>::iterator mit;
for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
ReadNC::VarData& vd = (*mit).second;
vd.entLoc = ReadNC::ENTLOCSET;
if (std::find(vd.varDims.begin(), vd.varDims.end(), tDim) != vd.varDims.end()) {
if ((std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) != vd.varDims.end()) &&
(std::find(vd.varDims.begin(), vd.varDims.end(), jCDim) != vd.varDims.end()))
vd.entLoc = ReadNC::ENTLOCFACE;
}
vd.numLev = 1;
if (std::find(vd.varDims.begin(), vd.varDims.end(), levDim) != vd.varDims.end())
vd.numLev = nLevels;
}
// For Eul models, slon and slat are "virtual" dimensions (not defined in the file header)
std::vector<std::string> ijdimNames(4);
ijdimNames[0] = "__slon";
ijdimNames[1] = "__slat";
ijdimNames[2] = "__lon";
ijdimNames[3] = "__lat";
std::string tag_name;
Tag tagh;
// __<dim_name>_LOC_MINMAX (for virtual slon, virtual slat, lon and lat)
for (unsigned int i = 0; i != ijdimNames.size(); i++) {
std::vector<int> val(2, 0);
if (ijdimNames[i] == "__slon") {
val[0] = lDims[0];
val[1] = lDims[3];
}
else if (ijdimNames[i] == "__slat") {
val[0] = lDims[1];
val[1] = lDims[4];
}
else if (ijdimNames[i] == "__lon") {
val[0] = lCDims[0];
val[1] = lCDims[3];
}
else if (ijdimNames[i] == "__lat") {
val[0] = lCDims[1];
val[1] = lCDims[4];
}
std::stringstream ss_tag_name;
ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
tag_name = ss_tag_name.str();
rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
ERRORR(rval, "Trouble creating __<dim_name>_LOC_MINMAX tag.");
rval = mbImpl->tag_set_data(tagh, &_fileSet, 1, &val[0]);
ERRORR(rval, "Trouble setting data for __<dim_name>_LOC_MINMAX tag.");
if (MB_SUCCESS == rval)
dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
}
// __<dim_name>_LOC_VALS (for virtual slon, virtual slat, lon and lat)
for (unsigned int i = 0; i != ijdimNames.size(); i++) {
void* val = NULL;
int val_len = 0;
if (ijdimNames[i] == "__slon") {
val = &ilVals[0];
val_len = ilVals.size();
}
else if (ijdimNames[i] == "__slat") {
val = &jlVals[0];
val_len = jlVals.size();
}
else if (ijdimNames[i] == "__lon") {
val = &ilCVals[0];
val_len = ilCVals.size();
}
else if (ijdimNames[i] == "__lat") {
val = &jlCVals[0];
val_len = jlCVals.size();
}
DataType data_type;
// Assume all has same data type as lon
switch (varInfo["lon"].varDataType) {
case NC_BYTE:
case NC_CHAR:
case NC_DOUBLE:
data_type = MB_TYPE_DOUBLE;
break;
case NC_FLOAT:
data_type = MB_TYPE_DOUBLE;
break;
case NC_INT:
data_type = MB_TYPE_INTEGER;
break;
case NC_SHORT:
default:
std::cerr << "Unrecognized data type for tag " << tag_name << std::endl;
ERRORR(MB_FAILURE, "Unrecognized data type");
break;
}
std::stringstream ss_tag_name;
ss_tag_name << ijdimNames[i] << "_LOC_VALS";
tag_name = ss_tag_name.str();
rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, data_type, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
ERRORR(rval, "Trouble creating __<dim_name>_LOC_VALS tag.");
rval = mbImpl->tag_set_by_ptr(tagh, &_fileSet, 1, &val, &val_len);
ERRORR(rval, "Trouble setting data for __<dim_name>_LOC_VALS tag.");
if (MB_SUCCESS == rval)
dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
}
// __<dim_name>_GLOBAL_MINMAX (for virtual slon, virtual slat, lon and lat)
for (unsigned int i = 0; i != ijdimNames.size(); i++) {
std::vector<int> val(2, 0);
if (ijdimNames[i] == "__slon") {
val[0] = gDims[0];
val[1] = gDims[3];
}
else if (ijdimNames[i] == "__slat") {
val[0] = gDims[1];
val[1] = gDims[4];
}
else if (ijdimNames[i] == "__lon") {
val[0] = gCDims[0];
val[1] = gCDims[3];
}
else if (ijdimNames[i] == "__lat") {
val[0] = gCDims[1];
val[1] = gCDims[4];
}
std::stringstream ss_tag_name;
ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
tag_name = ss_tag_name.str();
rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
ERRORR(rval, "Trouble creating __<dim_name>_GLOBAL_MINMAX tag.");
rval = mbImpl->tag_set_data(tagh, &_fileSet, 1, &val[0]);
ERRORR(rval, "Trouble setting data for __<dim_name>_GLOBAL_MINMAX tag.");
if (MB_SUCCESS == rval)
dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
}
// Hack: create dummy variables, if needed, for dimensions with no corresponding coordinate variables
rval = create_dummy_variables();
ERRORR(rval, "Failed to create dummy variables.");
return MB_SUCCESS;
}
1.7.6.1