Commit 97608a5f authored by David Thompson's avatar David Thompson Committed by Kitware Robot

Merge topic 'clm-pfb-files'

16ac8d26 Read ParFlow-formatted CLM state (`.C.pfb` files).
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: Sebastien Jourdain's avatarSebastien Jourdain <sebastien.jourdain@kitware.com>
Merge-request: !3237
parents 0f5c6020 16ac8d26
Pipeline #136116 running with stage
# ParFlow simulation support
* ParaView can now read ParFlow's simulation output
(files named `.pfb` and `.C.pfb`).
* Preliminary support for reading multiple ParFlow
files onto a single pair of meshes (surface and
subsurface) is also provided; it expects JSON
files with a `.pfmetadata` extension.
......@@ -21,6 +21,35 @@
>
<FileListDomain name="files"/>
</StringVectorProperty>
<IntVectorProperty
name="IsCLMFile"
label="Does file hold CLM data?"
animateable="1"
command="SetIsCLMFile"
default_values="-1"
number_of_elements="1"
>
<EnumerationDomain name="enum">
<Entry text="Infer from filename" value="-1"/>
<Entry text="No" value="0"/>
<Entry text="Yes" value="1"/>
</EnumerationDomain>
</IntVectorProperty>
<IntVectorProperty
name="CLMIrrType"
label="CLM irrigation type"
animateable="1"
command="SetCLMIrrType"
default_values="0"
number_of_elements="1"
>
<EnumerationDomain name="enum">
<Entry text="None" value="0"/>
<Entry text="Spray" value="1"/>
<Entry text="Drip" value="2"/>
<Entry text="Instant" value="3"/>
</EnumerationDomain>
</IntVectorProperty>
<Hints>
<ReaderFactory
extensions="pfb"
......@@ -45,6 +74,10 @@
proxygroup="internal_sources"
proxyname="ParFlowBinaryReader">
</Proxy>
<ExposedProperties>
<Property name="IsCLMFile"/>
<Property name="CLMIrrType"/>
</ExposedProperties>
</SubProxy>
<StringVectorProperty
......
......@@ -17,15 +17,39 @@
static constexpr std::streamoff headerSize = 6 * sizeof(double) + 4 * sizeof(int);
static constexpr std::streamoff subgridHeaderSize = 9 * sizeof(int);
static constexpr std::streamoff pfbEntrySize = sizeof(double);
static constexpr std::streamoff clmEntrySize =
11 * sizeof(double); // TODO: This is just a starting point.
static constexpr int clmBaseComponents = 11;
static constexpr std::streamoff clmEntrySize = clmBaseComponents * sizeof(double);
static const char* clmBaseComponentNames[clmBaseComponents] = { "eflx_lh_tot", "eflx_lwrad_out",
"eflx_sh_tot", "eflx_soil_grnd", "qflx_evap_tot", "qflx_evap_grnd", "qflx_evap_soi",
"qflx_evap_veg", "qflx_infl", "swe_out", "t_grnd" };
static std::streamoff computeEntrySize(bool isCLM, int nz)
{
std::streamoff sz;
if (!isCLM)
{
sz = pfbEntrySize;
}
else
{
if (nz < clmBaseComponents)
{
std::cerr << "Invalid CLM file: k axis must have at least " << clmBaseComponents
<< " entries but has " << nz << "\n";
}
sz = nz * sizeof(double);
}
return sz;
}
vtkStandardNewMacro(vtkParFlowReader);
vtkParFlowReader::vtkParFlowReader()
: FileName(nullptr)
, IsCLMFile(0)
, CLMIrrType(1)
, IsCLMFile(-1)
, CLMIrrType(0)
, NZ(0)
, InferredAsCLM(-1)
{
this->SetNumberOfInputPorts(0);
}
......@@ -39,7 +63,9 @@ void vtkParFlowReader::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os, indent);
os << indent << "FileName: " << (this->FileName ? this->FileName : "(null)") << "\n";
os << indent << "IsCLMFile: " << (this->IsCLMFile ? "true" : "false") << "\n";
os << indent
<< "IsCLMFile: " << (this->IsCLMFile > 0 ? "true" : this->IsCLMFile < 0 ? "infer" : "false")
<< "\n";
os << indent << "CLMIrrType: " << this->CLMIrrType << "\n";
os << indent << "IJKDivs:\n";
vtkIndent i2 = indent.GetNextIndent();
......@@ -93,11 +119,16 @@ int vtkParFlowReader::RequestData(vtkInformation* vtkNotUsed(request),
fnparts.push_back(token);
}
std::string arrayName("data");
this->InferredAsCLM = 0;
if (!fnparts.empty())
{
fnparts.pop_back(); // Drop the ".pfb".
if (!fnparts.empty())
{
if (fnparts.back() == "C" || fnparts.back() == "c")
{
this->InferredAsCLM = 1;
}
if (is_number(fnparts.back()))
{
fnparts.pop_back(); // Drop the timestep.
......@@ -108,6 +139,10 @@ int vtkParFlowReader::RequestData(vtkInformation* vtkNotUsed(request),
}
}
}
if (this->IsCLMFile >= 0)
{
this->InferredAsCLM = this->IsCLMFile;
}
// When run in parallel, we choose a range of blocks
// to load from those available.
......@@ -137,6 +172,8 @@ int vtkParFlowReader::RequestData(vtkInformation* vtkNotUsed(request),
vtkByteSwap::SwapBERange(dx.GetData(), 3);
vtkByteSwap::SwapBE(&numSubGrids);
this->NZ = nn[2]; // For computing size of CLM state.
std::streamoff measuredHeaderSize = pfb.tellg();
if (measuredHeaderSize != headerSize)
{
......@@ -168,6 +205,10 @@ int vtkParFlowReader::RequestData(vtkInformation* vtkNotUsed(request),
this->ReadBlock(pfb, output, xx, dx, arrayName, ni);
}
// Prevent accidents; don't preserve across calls to RequestData:
this->NZ = 0;
this->InferredAsCLM = -1;
return 1;
}
......@@ -265,7 +306,7 @@ std::streamoff vtkParFlowReader::GetBlockOffset(int blockId) const
std::streamoff vtkParFlowReader::GetBlockOffset(const vtkVector3i& blockIJK) const
{
std::streamoff offset;
std::streamoff entrySize = this->IsCLMFile ? clmEntrySize : pfbEntrySize;
std::streamoff entrySize = computeEntrySize(this->InferredAsCLM, this->NZ);
int ni = static_cast<int>(this->IJKDivs[0].size()) - 1;
int nj = static_cast<int>(this->IJKDivs[1].size()) - 1;
int blockId = blockIJK[0] + ni * (blockIJK[1] + nj * blockIJK[2]);
......@@ -291,7 +332,7 @@ std::streamoff vtkParFlowReader::GetBlockOffset(const vtkVector3i& blockIJK) con
std::streamoff vtkParFlowReader::GetEndOffset() const
{
std::streamoff offset;
std::streamoff entrySize = this->IsCLMFile ? clmEntrySize : pfbEntrySize;
std::streamoff entrySize = computeEntrySize(this->InferredAsCLM, this->NZ);
int ni = static_cast<int>(this->IJKDivs[0].size()) - 1;
int nj = static_cast<int>(this->IJKDivs[1].size()) - 1;
int nk = static_cast<int>(this->IJKDivs[2].size()) - 1;
......@@ -320,18 +361,82 @@ void vtkParFlowReader::ReadBlock(std::ifstream& pfb, vtkMultiBlockDataSet* outpu
if (this->ReadSubgridHeader(pfb, si, sn, sr))
{
vtkNew<vtkImageData> image;
vtkNew<vtkDoubleArray> field;
image->SetOrigin(origin.GetData());
image->SetSpacing(spacing.GetData());
image->SetExtent(si[0], si[0] + sn[0], si[1], si[1] + sn[1], si[2], si[2] + sn[2]);
field->SetNumberOfTuples(image->GetNumberOfCells());
field->SetName(arrayName.c_str());
image->GetCellData()->SetScalars(field);
vtkIdType numValues = sn[0] * sn[1] * sn[2];
pfb.read(reinterpret_cast<char*>(field->GetVoidPointer(0)), sizeof(double) * numValues);
vtkByteSwap::SwapBERange(reinterpret_cast<double*>(field->GetVoidPointer(0)), numValues);
if (this->InferredAsCLM)
{
// The CLM files have the full simulation extent listed but only
// provide data on the top 2-d surface:
image->SetExtent(si[0], si[0] + sn[0], si[1], si[1] + sn[1], si[2], si[2]);
const int numCLMVars = sn[2] - si[2];
int numComponents = 0;
for (int cc = 0; cc < clmBaseComponents && cc < numCLMVars; ++cc, ++numComponents)
{
vtkNew<vtkDoubleArray> field;
field->SetName(clmBaseComponentNames[cc]);
this->ReadBlockIntoArray(pfb, image, field);
}
switch (this->CLMIrrType)
{
case 1:
{
vtkNew<vtkDoubleArray> field;
field->SetName("qflx_qirr");
this->ReadBlockIntoArray(pfb, image, field);
++numComponents;
}
break;
case 3:
{
vtkNew<vtkDoubleArray> field;
field->SetName("qflx_qirr_inst");
this->ReadBlockIntoArray(pfb, image, field);
++numComponents;
}
break;
default:
break;
}
for (int cz = 0; numComponents < numCLMVars; ++cz, ++numComponents)
{
vtkNew<vtkDoubleArray> field;
std::ostringstream name;
name << "tsoil_" << cz;
field->SetName(name.str().c_str());
this->ReadBlockIntoArray(pfb, image, field);
}
}
else
{
// Read a single PFB state variable:
vtkNew<vtkDoubleArray> field;
image->SetExtent(si[0], si[0] + sn[0], si[1], si[1] + sn[1], si[2], si[2] + sn[2]);
field->SetName(arrayName.c_str());
this->ReadBlockIntoArray(pfb, image, field);
}
output->SetBlock(blockId, image);
}
}
void vtkParFlowReader::ReadBlockIntoArray(
std::ifstream& file, vtkImageData* img, vtkDoubleArray* arr)
{
arr->SetNumberOfTuples(img->GetNumberOfCells());
auto cellData = img->GetCellData();
// Calling cellData->SetScalars(arr) multiple times removes
// previous arrays set as scalars, so be careful not to:
cellData->AddArray(arr);
if (!cellData->GetScalars())
{
cellData->SetScalars(arr);
}
vtkIdType numValues = arr->GetNumberOfTuples() * arr->GetNumberOfComponents();
file.read(reinterpret_cast<char*>(arr->GetVoidPointer(0)), sizeof(double) * numValues);
vtkByteSwap::SwapBERange(reinterpret_cast<double*>(arr->GetVoidPointer(0)), numValues);
// std::cout << arr->GetNumberOfComponents() << " " << numValues << "Read to byte " << pfb.tellg()
// << "\n";
}
......@@ -10,6 +10,8 @@
#include <fstream>
#include <vector>
class vtkDoubleArray;
class vtkImageData;
class vtkMultiBlockDataSet;
/**\brief Read ParFlow simulation output.
......@@ -21,6 +23,13 @@ class vtkMultiBlockDataSet;
* If there are fewer blocks than ranks, some processes will do no work.
* You may use "pftools dist" to repartition the data into a different
* number of blocks (known in ParFlow as subgrids).
*
* You may indicate files contain CLM data or a single PFB state variable.
* CLM files have multiple values per cell while PFB files hold only one.
* The per-cell values are stored as i-j image stacks (i.e., each subgrid
* stores a sequence of 2-d images, one per CLM state variable); the k-index
* extent of the PFB file corresponds to the number of CLM state variables
* per cell.
*/
class VTKPARFLOWIO_EXPORT vtkParFlowReader : public vtkMultiBlockDataSetAlgorithm
{
......@@ -35,17 +44,21 @@ public:
/// Set/get whether the file is a CLM output (i.e., ".c.pfb").
///
/// The default is false.
/// CLM files have multiple values per cell while PFB files have only one.
/// This variable accepts 3 values:
/// + 0 indicates the file is not a CLM file
/// + 1 indicates the file is a CLM file
/// + -1 indicates the filetype should be inferred from the filename.
///
/// The default is -1.
vtkGetMacro(IsCLMFile, int);
vtkSetMacro(IsCLMFile, int);
vtkBooleanMacro(IsCLMFile, int);
/// Set/get the CLM irradiance type.
/// Set/get the CLM irrigation type.
///
/// The CLM irradiance types are enumerated by ParFlow.
/// At least 1 and 3 are valid values.
/// The default is 1.
/// The CLM irrigation types are enumerated by ParFlow:
/// 0=None, 1=Spray, 2=Drip, 3=Instant.
/// The default is 0.
vtkGetMacro(CLMIrrType, int);
vtkSetMacro(CLMIrrType, int);
......@@ -91,12 +104,17 @@ protected:
static void GetBlockExtent(const vtkVector3i& wholeExtentIn, const vtkVector3i& numberOfBlocksIn,
const vtkVector3i& blockIJKIn, vtkVector3i& blockExtentMinOut, vtkVector3i& blockExtentMaxOut);
static void ReadBlockIntoArray(std::ifstream& file, vtkImageData* img, vtkDoubleArray* arr);
/// The filename, which must be a valid path before RequestData is called.
char* FileName;
int IsCLMFile;
int CLMIrrType;
/// IJKDiv only valid inside RequestData; used to compute subgrid offsets.
/// IJKDivs, NZ, and InferredAsCLM are only valid inside RequestData; used to compute subgrid
/// offsets.
std::vector<int> IJKDivs[3];
int NZ;
int InferredAsCLM;
};
#endif // vtkParflowReader_h
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment