6

对于盖驱动腔 (CFD) 的格子玻尔兹曼模拟,我将立方域分解为(也是立方的)8 个子域,这些子域由 8 个等级独立计算。每个 MPI 等级都会为每个时间步生成一个 VTK 文件,由于我使用的是 ParaView,我想将整个事物可视化为一个立方体。更具体地说明我想要实现的目标:

  • 我有一个长度为 8 的立方体(每个方向的元素数)=> 8x8x8 = 512 个元素。
  • 每个维度分布到 2 个等级,即每个等级处理 4x4x4 = 64 个元素。
  • 每个等级都将其结果写入lbm_out_<rank>.<timestep>.vts以 VTKStructuredGrid格式调用的文件。
  • 我想生成一个.pvts文件,该文件收集*.vts文件并将包含子域的文件组合到 ParaView 可以视为整个域的单个文件中。

不幸的是,我遇到了很多问题,因为我觉得 ParaView 和 VTK 的文档记录极差,而且来自 ParaView 的错误消息完全没用。

我有以下*.pvts文件,其中包括一个幽灵层和:

<?xml version="1.0"?>
<VTKFile type="PStructuredGrid" version="0.1" byte_order="LittleEndian">
    <PStructuredGrid WholeExtent="0 7 0 7 0 7 " GhostLevel="1">
        <PPoints>
            <PDataArray NumberOfComponents="3" type="Float32" />
        </PPoints>
        <Piece Extent="0 4 0 4 0 4" Source="lbm_out_0.0.vts"/>
        <Piece Extent="3 7 0 4 0 4" Source="lbm_out_1.0.vts"/>
        <Piece Extent="0 4 3 7 0 4" Source="lbm_out_2.0.vts"/>
        <Piece Extent="3 7 3 7 0 4" Source="lbm_out_3.0.vts"/>
        <Piece Extent="0 4 0 4 3 7" Source="lbm_out_4.0.vts"/>
        <Piece Extent="3 7 0 4 3 7" Source="lbm_out_5.0.vts"/>
        <Piece Extent="0 4 3 7 3 7" Source="lbm_out_6.0.vts"/>
        <Piece Extent="3 7 3 7 3 7" Source="lbm_out_7.0.vts"/>
    </PStructuredGrid>
</VTKFile>

有了那个我觉得应该可以正常工作的文件(请注意,还没有参数,只是简单的几何信息),我的域范围完全搞砸了,尽管每个*.vts文件都可以自己正常工作。我附上了一个屏幕截图,以使事情更清楚: 请注意界限是如何完全超出范围的。

可能是什么问题?是否可以将旧版 VTK 文件用于此任务?我可能做错了什么吗?我真的不知道如何完成这项任务,而且我在 google 中找到的资源非常有限。谢谢你。

4

1 回答 1

4

不幸的是,没有vtkXMLPStructuredGridWriter的示例(示例中未使用 VTK 类)。所以我决定编写最简单的代码来为结构化网格生成 *.vts 和 .pvts 文件,这与您正在寻找的情况非常相似。

以下代码使用 MPI 和 VTK 编写并行结构化网格文件。在此示例中,我们有两个进程创建自己的 .vts 文件,vtkXMLPStructuredGridWriter 类写入 .pvts 文件:

// MPI Library
#include <mpi.h>

//VTK Library
#include <vtkXMLPStructuredGridWriter.h>
#include <vtkStructuredGrid.h>
#include <vtkSmartPointer.h>
#include <vtkFloatArray.h>
#include <vtkCellData.h>
#include <vtkMPIController.h>
#include <vtkProgrammableFilter.h>
#include <vtkInformation.h>

struct Args {
  vtkProgrammableFilter* pf;
  int local_extent[6];
};

// function to operate on the point attribute data
void execute (void* arg) {
  Args* args = reinterpret_cast<Args*>(arg);
  auto info = args->pf->GetOutputInformation(0);
  auto output_tmp = args->pf->GetOutput();
  auto input_tmp  = args->pf->GetInput();
  vtkStructuredGrid* output = dynamic_cast<vtkStructuredGrid*>(output_tmp);
  vtkStructuredGrid* input  = dynamic_cast<vtkStructuredGrid*>(input_tmp);
  output->ShallowCopy(input);
  output->SetExtent(args->local_extent);
}

int main (int argc, char *argv[]) {
  MPI_Init (&argc, &argv);
  int myrank;
  MPI_Comm_rank (MPI_COMM_WORLD, &myrank);

  int lx {10}, ly{10}, lz{10};        //local dimensions of the process's grid
  int dims[3] = {lx+1, ly+1, lz+1};
  int global_extent[6] = {0, 2*lx, 0, ly, 0, lz};
  int local_extent[6] = {myrank*lx, (myrank+1)*lx, 0, ly, 0, lz};

  // Create and Initialize vtkMPIController
  auto contr = vtkSmartPointer<vtkMPIController>::New();
  contr->Initialize(&argc, &argv, 1);
  int nranks = contr->GetNumberOfProcesses();
  int rank   = contr->GetLocalProcessId();

  // Create grid points, allocate memory and Insert them
  auto points = vtkSmartPointer<vtkPoints>::New();
  points->Allocate(dims[0]*dims[1]*dims[2]);
  for (int k=0; k<dims[2]; ++k)
    for (int j=0; j<dims[1]; ++j)
      for (int i=0; i<dims[0]; ++i)
        points->InsertPoint(i + j*dims[0] + k*dims[0]*dims[1],
                            i+rank*(dims[0]-1), j, k);

  // Create a density field. Note that the number of cells is always less than
  // number of grid points by an amount of one so we use dims[i]-1
  auto density = vtkSmartPointer<vtkFloatArray>::New();
  density->SetNumberOfComponents(1);
  density->SetNumberOfTuples((dims[0]-1)*(dims[1]-1)*(dims[2]-1));
  density->SetName ("density");
  int index;
  for (int k=0; k<lz; ++k)
    for (int j=0; j<ly; ++j)
      for (int i=0; i<lx; ++i) {
        index = i + j*(dims[0]-1) + k*(dims[0]-1)*(dims[1]-1);
        density->SetValue(index, i+j+k);
      }

  // Create a vtkProgrammableFilter
  auto pf = vtkSmartPointer<vtkProgrammableFilter>::New();

  // Initialize an instance of Args
  Args args;
  args.pf = pf;
  for(int i=0; i<6; ++i) args.local_extent[i] = local_extent[i];

  pf->SetExecuteMethod(execute, &args);

  // Create a structured grid and assign point data and cell data to it
  auto structuredGrid = vtkSmartPointer<vtkStructuredGrid>::New();
  structuredGrid->SetExtent(global_extent);
  pf->SetInputData(structuredGrid);
  structuredGrid->SetPoints(points);
  structuredGrid->GetCellData()->AddArray(density);

  // Create the parallel writer and call some functions
  auto parallel_writer = vtkSmartPointer<vtkXMLPStructuredGridWriter>::New();
  parallel_writer->SetInputConnection(pf->GetOutputPort());
  parallel_writer->SetController(contr);
  parallel_writer->SetFileName("test.pvts");
  parallel_writer->SetNumberOfPieces(nranks);
  parallel_writer->SetStartPiece(rank);
  parallel_writer->SetEndPiece(rank);
  parallel_writer->SetDataModeToBinary();
  parallel_writer->Update();
  parallel_writer->Write();

  contr->Finalize();

  // WARNING: it seems that MPI_Finalize is not necessary since we are using
  // Finalize() function from vtkMPIController class. Uncomment the following
  // line and see what happens.
//   MPI_Finalize ();
  return 0;
}

此代码仅将一些数据(在本例中为标量的密度)写入文件而不渲染它。对于可视化,您将需要像 Paraview 这样的软件。要运行此代码,您可以使用此 CMake 文件:

cmake_minimum_required(VERSION 2.8)

PROJECT(PXMLStructuredGridWriter)
add_executable(PXMLStructuredGridWriter parallel_vtk.cpp)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})
target_link_libraries(PXMLStructuredGridWriter ${VTK_LIBRARIES})

find_package(MPI REQUIRED)
include_directories(${MPI_INCLUDE_PATH})
target_link_libraries(PXMLStructuredGridWriter ${MPI_LIBRARIES})

最后,您会在与可执行文件相同的目录中找到这样的 xml 文件:

<?xml version="1.0"?>
<VTKFile type="PStructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
  <PStructuredGrid WholeExtent="0 20 0 10 0 10" GhostLevel="2">
    <PCellData>
      <PDataArray type="Float32" Name="density"/>
    </PCellData>
    <PPoints>
      <PDataArray type="Float32" Name="Points" NumberOfComponents="3"/>
    </PPoints>
    <Piece Extent="0 10 0 10 0 10" Source="test_0.vts"/>
    <Piece Extent="10 20 0 10 0 10" Source="test_1.vts"/>
  </PStructuredGrid>
</VTKFile>
于 2018-04-25T10:26:31.590 回答