moab
|
Child helper class for Finite Volume grid (CAM_FV) More...
#include <NCHelperFV.hpp>
Public Member Functions | |
NCHelperFV (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 Finite Volume grid (CAM_FV)
Definition at line 17 of file NCHelperFV.hpp.
moab::NCHelperFV::NCHelperFV | ( | ReadNC * | readNC, |
int | fileId, | ||
const FileOptions & | opts, | ||
EntityHandle | fileSet | ||
) | [inline] |
Definition at line 20 of file NCHelperFV.hpp.
: ScdNCHelper(readNC, fileId, opts, fileSet) {}
bool moab::NCHelperFV::can_read_file | ( | ReadNC * | readNC, |
int | fileId | ||
) | [static] |
Definition at line 13 of file NCHelperFV.cpp.
{ std::vector<std::string>& dimNames = readNC->dimNames; // 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("lon")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end())) { // 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::NCHelperFV::get_mesh_type_name | ( | ) | [inline, private, virtual] |
ErrorCode moab::NCHelperFV::init_mesh_vals | ( | ) | [private, virtual] |
Interfaces to be implemented in child classes.
Implements moab::NCHelper.
Definition at line 45 of file NCHelperFV.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 i/j dimensions // First i std::vector<std::string>::iterator vit; unsigned int idx; if ((vit = std::find(dimNames.begin(), dimNames.end(), "slon")) != dimNames.end()) idx = vit - dimNames.begin(); else { ERRORR(MB_FAILURE, "Couldn't find 'slon' variable."); } iDim = idx; gDims[0] = 0; gDims[3] = dimLens[idx] - 1; // Then j if ((vit = std::find(dimNames.begin(), dimNames.end(), "slat")) != dimNames.end()) idx = vit - dimNames.begin(); else { ERRORR(MB_FAILURE, "Couldn't find 'slat' variable."); } jDim = idx; gDims[1] = 0; gDims[4] = dimLens[idx] - 1 + 2; // Add 2 for the pole points // Look for names of center i/j dimensions // First i if ((vit = std::find(dimNames.begin(), dimNames.end(), "lon")) != dimNames.end()) idx = vit - dimNames.begin(); else { ERRORR(MB_FAILURE, "Couldn't find 'lon' variable."); } 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; if (globallyPeriodic[0]) assert("Number of vertices and edges should be same" && gDims[3] == gCDims[3]); else assert("Number of vertices should equal to number of edges plus one" && gDims[3] == gCDims[3] + 1); // 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 FV models, will always be non-periodic in j assert(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]; // Get number of levels 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 FV 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("slon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) { // Last column if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3]) { assert(lDims[3] == gDims[3] + 1); std::vector<double> dummyVar(lDims[3] - lDims[0]); rval = read_coordinate("slon", lDims[0], lDims[3] - 1, dummyVar); double dif = dummyVar[1] - dummyVar[0]; std::size_t i; for (i = 0; i != dummyVar.size(); i++) ilVals[i] = dummyVar[i]; ilVals[i] = ilVals[i - 1] + dif; } else { rval = read_coordinate("slon", lDims[0], lDims[3], ilVals); ERRORR(rval, "Trouble reading 'slon' variable."); } } else { ERRORR(MB_FAILURE, "Couldn't find 'slon' variable."); } } if (-1 != lDims[1]) { if ((vmit = varInfo.find("slat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) { if (!isParallel || ((gDims[4] - gDims[1]) == (lDims[4] - lDims[1]))) { std::vector<double> dummyVar(lDims[4] - lDims[1] - 1); rval = read_coordinate("slat", lDims[1], lDims[4] - 2, dummyVar); ERRORR(rval, "Trouble reading 'slat' variable."); // Copy the correct piece jlVals[0] = -90.0; std::size_t i = 0; for (i = 1; i != dummyVar.size() + 1; i++) jlVals[i] = dummyVar[i - 1]; jlVals[i] = 90.0; // Using value of i after loop exits. } else { // If this is the first row // Need to read one less then available and read it into a dummy var if (lDims[1] == gDims[1]) { std::vector<double> dummyVar(lDims[4] - lDims[1]); rval = read_coordinate("slat", lDims[1], lDims[4] - 1, dummyVar); ERRORR(rval, "Trouble reading 'slat' variable."); // Copy the correct piece jlVals[0] = -90.0; for (int i = 1; i < lDims[4] + 1; i++) jlVals[i] = dummyVar[i - 1]; } // Or if it's the last row else if (lDims[4] == gDims[4]) { std::vector<double> dummyVar(lDims[4] - lDims[1]); rval = read_coordinate("slat", lDims[1] - 1, lDims[4] - 2, dummyVar); ERRORR(rval, "Trouble reading 'slat' variable."); // Copy the correct piece std::size_t i = 0; for (i = 0; i != dummyVar.size(); i++) jlVals[i] = dummyVar[i]; jlVals[i] = 90.0; // Using value of i after loop exits. } // It's in the middle else { rval = read_coordinate("slat", lDims[1] - 1, lDims[4] - 1, jlVals); ERRORR(rval, "Trouble reading 'slat' variable."); } } } else { ERRORR(MB_FAILURE, "Couldn't find 'slat' 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; else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) != vd.varDims.end())) vd.entLoc = ReadNC::ENTLOCNSEDGE; else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jCDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(), vd.varDims.end(), iDim) != vd.varDims.end())) vd.entLoc = ReadNC::ENTLOCEWEDGE; } vd.numLev = 1; if (std::find(vd.varDims.begin(), vd.varDims.end(), levDim) != vd.varDims.end()) vd.numLev = nLevels; } 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 slon, 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 slon, 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 slon, 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; }