Description
I am using custom solver to work. I use VTK to output VTS file, when three dimensions are the same, there is no errors. But when they are not the same, paraview read the file and shows weird wireframe preview, looks like several line overlapping. What is more, there is ripple-like disappearance of values along certain dimension, such as velocity.
I output the same fluid grid to csv, and then transform it into vts using pyevtk. The result is correct. No error in paraview.
So there must be some different between my output class code and others, such as pyevtk. But I don’t know what is it. Any idea? Thank you!
Figures
For example, it is resample vorticity to image at time 0, looking from X axis. Using pvNVIDIAIndex plugin.
After steps of iteration, the fluid becomes turbulent. But the disappearance of values remains.
So it is not my NS solver’s fault.
It is weird wireframe preview.
And I have tried to reduce my dimension to 20, 10, 10. The I keep showing my vts output, toggling the pyevtk output, recording the screen. You can see the right output wireframe by pyevtk is the bounding box of my error vtk wireframe.
Codes
My part of code related about vts output is:
#pragma once
#include "fluid_grid/fluid_grid_3d.h"
#include <cassert>
#include <vtkDoubleArray.h>
#include <vtkErrorCode.h>
#include <vtkFloatArray.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkSmartPointer.h>
#include <vtkStructuredGrid.h>
#include <vtkXMLStructuredGridWriter.h>
class VTSWriter
{
private:
vtkSmartPointer<vtkPoints> points;
public:
VTSWriter(int Nx, int Ny, int Nz, double h, double offset_x = 0.0, double offset_y = 0.0, double offset_z = 0.0)
{
assert(Nx > 0);
assert(Ny > 0);
assert(Nz > 0);
assert(h > 0.0);
points = vtkSmartPointer<vtkPoints>::New();
for (int i = 0; i < Nx; i++)
{
for (int j = 0; j < Ny; j++)
{
for (int k = 0; k < Nz; k++)
{
double x = offset_x + i * h + 0.5 * h;
double y = offset_y + j * h + 0.5 * h;
double z = offset_z + k * h + 0.5 * h;
points->InsertNextPoint(x, y, z);
}
}
}
}
void write(FluidGrid3D& grid, const std::string& filename);
};
#include "vts_writer.h"
#include "file_monitor.h"
void VTSWriter::write(FluidGrid3D& grid, const std::string& filename)
{
EXPOSE_GRID3D_PARA(grid);
auto& u = grid.u;
auto& v = grid.v;
auto& w = grid.w;
auto& p = grid.p;
vtkSmartPointer<vtkStructuredGrid> vtk_structured_grid = vtkSmartPointer<vtkStructuredGrid>::New();
vtkSmartPointer<vtkDoubleArray> velocity = vtkSmartPointer<vtkDoubleArray>::New();
velocity->SetName("Velocity");
velocity->SetNumberOfComponents(3);
velocity->SetNumberOfTuples(Nx * Ny * Nz);
vtkSmartPointer<vtkDoubleArray> pressure = vtkSmartPointer<vtkDoubleArray>::New();
pressure->SetName("Pressure");
pressure->SetNumberOfComponents(1);
pressure->SetNumberOfTuples(Nx * Ny * Nz);
int idx = 0;
for (int i = 0; i < Nx; i++)
{
for (int j = 0; j < Ny; j++)
{
for (int k = 0; k < Nz; k++)
{
double uc = (u(i + usi - 1, j + usj, k + usk) + u(i + usi, j + usj, k + usk)) / 2.0;
double vc = (v(i + vsi, j + vsj - 1, k + vsk) + v(i + vsi, j + vsj, k + vsk)) / 2.0;
double wc = (w(i + wsi, j + wsj, k + wsk - 1) + w(i + wsi, j + wsj, k + wsk)) / 2.0;
double pc = p(i, j, k);
velocity->SetTuple3(idx, uc, vc, wc);
pressure->SetTuple1(idx, pc);
idx++;
}
}
}
vtk_structured_grid->SetDimensions(Nx, Ny, Nz);
vtk_structured_grid->SetPoints(points);
vtk_structured_grid->GetPointData()->AddArray(velocity);
vtk_structured_grid->GetPointData()->AddArray(pressure);
std::atomic<bool> stop_flag(false);
std::streampos expected_size = (std::size_t)Nx * Ny * Nz * (9 * 4 + 8);
std::thread monitor_thread(monitor_file_size, filename + ".vts", expected_size, std::ref(stop_flag));
vtkSmartPointer<vtkXMLStructuredGridWriter> writer = vtkSmartPointer<vtkXMLStructuredGridWriter>::New();
writer->SetDataModeToAppended();
writer->SetCompressorTypeToNone();
writer->EncodeAppendedDataOff();
writer->SetFileName((filename + ".vts").c_str());
writer->SetInputData(vtk_structured_grid);
writer->Write();
stop_flag = true;
monitor_thread.join();
if (writer->GetErrorCode() != vtkErrorCode::NoError)
{
std::cout << "Writing to file failed: " << filename + ".vts" << std::endl;
}
}
And my csv output code is
void fluidgrid_to_csv(FluidGrid3D& grid, const std::string& filename, double offset_x, double offset_y, double offset_z)
{
std::ofstream outfile(filename + ".csv");
if (!outfile.is_open())
{
std::cerr << "Failed to open file: " << filename + ".csv" << std::endl;
return;
}
std::cout << "Writing csv: " << filename + ".csv" << std::endl;
outfile << "x,y,z,u,v,w,p" << std::endl;
EXPOSE_GRID3D_PARA(grid);
auto& u = grid.u;
auto& v = grid.v;
auto& w = grid.w;
auto& p = grid.p;
std::atomic<bool> stop_flag(false);
std::streampos expected_size = (std::size_t)Nx * Ny * Nz * 70;
std::thread monitor_thread(monitor_file_size, filename + ".csv", expected_size, std::ref(stop_flag));
int precision = 6;
outfile << std::fixed << std::setprecision(precision);
for (int i = 0; i < Nx; i++)
{
for (int j = 0; j < Ny; j++)
{
for (int k = 0; k < Nz; k++)
{
// x, y, z
outfile << offset_x + i * h + 0.5 * h << "," << offset_y + j * h + 0.5 * h << ","
<< offset_z + k * h + 0.5 * h << ",";
double uc = (u(i + usi - 1, j + usj, k + usk) + u(i + usi, j + usj, k + usk)) / 2.0;
double vc = (v(i + vsi, j + vsj - 1, k + vsk) + v(i + vsi, j + vsj, k + vsk)) / 2.0;
double wc = (w(i + wsi, j + wsj, k + wsk - 1) + w(i + wsi, j + wsj, k + wsk)) / 2.0;
double pc = p(i, j, k);
// u, v, w
outfile << uc << ",";
outfile << vc << ",";
outfile << wc << ",";
// p
outfile << pc << std::endl;
}
}
}
outfile << std::endl;
outfile.close();
stop_flag = true;
monitor_thread.join();
if (outfile.fail())
{
std::cout << "Failed to write file: " << filename + ".csv" << std::endl;
}
}
You can see my output content are the same between csv writer and vts writer.
My test
I have writed a comparing test, to read two vts file and compare the value. Following code will output the first vts file’s velocity field, and then output the difference of velocity field between the two vts file.
#include "fluid_grid/fluid_grid_3d.h"
#include <iostream>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkSmartPointer.h>
#include <vtkStructuredGrid.h>
#include <vtkXMLStructuredGridReader.h>
void ReadVTS(FluidGrid3D& grid, const std::string& filename)
{
vtkSmartPointer<vtkXMLStructuredGridReader> reader = vtkSmartPointer<vtkXMLStructuredGridReader>::New();
reader->SetFileName(filename.c_str());
reader->Update();
vtkSmartPointer<vtkStructuredGrid> vtk_structured_grid = reader->GetOutput();
vtkSmartPointer<vtkDoubleArray> velocity =
vtkDoubleArray::SafeDownCast(vtk_structured_grid->GetPointData()->GetArray("Velocity"));
vtkSmartPointer<vtkDoubleArray> pressure =
vtkDoubleArray::SafeDownCast(vtk_structured_grid->GetPointData()->GetArray("Pressure"));
if (!velocity || !pressure)
{
std::cerr << "Error: Unable to find required data arrays in the VTS file." << std::endl;
return;
}
EXPOSE_GRID3D_PARA(grid);
for (int i = 0; i < Nx; i++)
{
for (int j = 0; j < Ny; j++)
{
for (int k = 0; k < Nz; k++)
{
int idx = i * Ny * Nz + j * Nz + k;
double uc, vc, wc, pc;
velocity->GetTuple(idx, &uc);
velocity->GetTuple(idx, &vc);
velocity->GetTuple(idx, &wc);
pc = pressure->GetTuple1(idx);
grid.u(i + usi - 1, j + usj, k + usk) = uc;
grid.v(i + vsi, j + vsj - 1, k + vsk) = vc;
grid.w(i + wsi, j + wsj, k + wsk - 1) = wc;
grid.p(i, j, k) = pc;
}
}
}
}
void Compare(FluidGrid3D& grid1, FluidGrid3D& grid2)
{
EXPOSE_GRID3D_PARA(grid1);
std::ofstream outfile("comp_result.csv");
if (!outfile.is_open())
{
std::cerr << "Failed to open file: " << "comp_result.csv" << std::endl;
return;
}
for (int i = 0; i < Nx; i++)
{
for (int j = 0; j < Ny; j++)
{
for (int k = 0; k < Nz; k++)
{
int idx = i * Ny * Nz + j * Nz + k;
outfile << grid1.u(i + usi - 1, j + usj, k + usk) << ",";
outfile << grid1.v(i + vsi, j + vsj - 1, k + vsk) << ",";
outfile << grid1.w(i + wsi, j + wsj, k + wsk - 1) << std::endl;
}
}
}
for (int i = 0; i < Nx; i++)
{
for (int j = 0; j < Ny; j++)
{
for (int k = 0; k < Nz; k++)
{
int idx = i * Ny * Nz + j * Nz + k;
outfile << grid1.u(i + usi - 1, j + usj, k + usk) - grid2.u(i + usi - 1, j + usj, k + usk) << ",";
outfile << grid1.v(i + vsi, j + vsj - 1, k + vsk) - grid2.v(i + vsi, j + vsj - 1, k + vsk) << ",";
outfile << grid1.w(i + wsi, j + wsj, k + wsk - 1) - grid2.w(i + wsi, j + wsj, k + wsk - 1) << std::endl;
}
}
}
outfile.close();
if (outfile.fail())
{
std::cout << "Writing to file failed: " << "comp_result.csv" << std::endl;
}
}
int main(int argc, char* argv[])
{
if (argc != 3)
{
std::cout << "Para number should be 3." << std::endl;
return 1;
}
constexpr double lx = 80.0;
constexpr int Nx = 20;
constexpr int Ny = 10;
constexpr int Nz = 10;
constexpr double H = lx / Nx;
FluidGrid3D fluid_grid1 = CREATE_GRID3D(Nx, Ny, Nz, H);
FluidGrid3D fluid_grid2 = CREATE_GRID3D(Nx, Ny, Nz, H);
ReadVTS(fluid_grid1, argv[1]);
ReadVTS(fluid_grid2, argv[1]);
Compare(fluid_grid1, fluid_grid2);
return 0;
}
I tested vts file created by VTK and by pyevtk, the result is no difference.
I also compare there position
#include <iostream>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkSmartPointer.h>
#include <vtkStructuredGrid.h>
#include <vtkXMLStructuredGridReader.h>
#include <fstream>
void ReadPositions(const std::string& filename, double*& x, double*& y, double*& z, int& num_points)
{
vtkSmartPointer<vtkXMLStructuredGridReader> reader = vtkSmartPointer<vtkXMLStructuredGridReader>::New();
reader->SetFileName(filename.c_str());
reader->Update();
vtkSmartPointer<vtkStructuredGrid> vtk_structured_grid = reader->GetOutput();
num_points = vtk_structured_grid->GetNumberOfPoints();
if (num_points <= 0)
{
std::cerr << "Error: No points found in the VTS file." << std::endl;
return;
}
x = new double[num_points];
y = new double[num_points];
z = new double[num_points];
for (int i = 0; i < num_points; i++)
{
x[i] = 0;
y[i] = 0;
z[i] = 0;
}
for (int i = 0; i < num_points; ++i)
{
double point[3];
vtk_structured_grid->GetPoint(i, point);
x[i] = point[0];
y[i] = point[1];
z[i] = point[2];
}
}
void ComparePositions(double* x1, double* y1, double* z1, double* x2, double* y2, double* z2, int num_points1, int num_points2)
{
if (num_points1 != num_points2)
{
std::cerr << "Error: The two VTS files have a different number of points." << std::endl;
return;
}
std::ofstream outfile("position_comparison.csv");
if (!outfile.is_open())
{
std::cerr << "Failed to open file: " << "position_comparison.csv" << std::endl;
return;
}
outfile << "Index,x1,y1,z1,x2,y2,z2,dx,dy,dz\n";
for (int i = 0; i < num_points1; ++i)
{
double dx = x1[i] - x2[i];
double dy = y1[i] - y2[i];
double dz = z1[i] - z2[i];
outfile << i << ","
<< x1[i] << "," << y1[i] << "," << z1[i] << ","
<< x2[i] << "," << y2[i] << "," << z2[i] << ","
<< dx << "," << dy << "," << dz << "\n";
}
outfile.close();
if (outfile.fail())
{
std::cerr << "Writing to file failed: " << "position_comparison.csv" << std::endl;
}
}
int main(int argc, char* argv[])
{
if (argc != 3)
{
std::cout << "Usage: " << argv[0] << " <filename1> <filename2>" << std::endl;
return 1;
}
double* x1 = nullptr;
double* y1 = nullptr;
double* z1 = nullptr;
int num_points1 = 0;
double* x2 = nullptr;
double* y2 = nullptr;
double* z2 = nullptr;
int num_points2 = 0;
ReadPositions(argv[1], x1, y1, z1, num_points1);
ReadPositions(argv[2], x2, y2, z2, num_points2);
if (x1 != nullptr && y1 != nullptr && z1 != nullptr && x2 != nullptr && y2 != nullptr && z2 != nullptr)
{
ComparePositions(x1, y1, z1, x2, y2, z2, num_points1, num_points2);
delete[] x1;
delete[] y1;
delete[] z1;
delete[] x2;
delete[] y2;
delete[] z2;
}
return 0;
}
But the result is strange, vts created by VTK has right position, but vts created by pyevtk has error position
Index,x1,y1,z1,x2,y2,z2,dx,dy,dz
0,2,2,2,0,5.30499e-315,5.30499e-315,2,2,2
1,2,2,6,5.30499e-315,5.31276e-315,5.30499e-315,2,2,6
2,2,2,10,5.30499e-315,5.31665e-315,5.30499e-315,2,2,10
3,2,2,14,5.30499e-315,5.31924e-315,5.30499e-315,2,2,14
4,2,2,18,5.30499e-315,5.32118e-315,5.30499e-315,2,2,18
5,2,2,22,5.30499e-315,5.32247e-315,5.30499e-315,2,2,22
x1,y1,z1 correspond to vts created by VTK, x2,y2,z2 correspond to vts created by pyevtk.
Minial Reproduction
#include <vtkNew.h>
#include <vtkPoints.h>
#include <vtkStructuredGrid.h>
#include <vtkXMLStructuredGridWriter.h>
int main(int, char*[])
{
// Create a grid.
vtkNew<vtkStructuredGrid> structuredGrid;
vtkNew<vtkPoints> points;
int Nx = 20;
int Ny = 10;
int Nz = 10;
double h = 0.1;
for (int i = 0; i < Nx; i++)
{
for (int j = 0; j < Ny; j++)
{
for (int k = 0; k < Nz; k++)
{
double x = 0 + i * h + 0.5 * h;
double y = 0 + j * h + 0.5 * h;
double z = 0 + k * h + 0.5 * h;
points->InsertNextPoint(x, y, z);
}
}
}
// Specify the dimensions of the grid.
structuredGrid->SetDimensions(Nx, Ny, Nz);
structuredGrid->SetPoints(points);
// Write file.
vtkNew<vtkXMLStructuredGridWriter> writer;
writer->SetFileName("output.vts");
writer->SetInputData(structuredGrid);
writer->Write();
return EXIT_SUCCESS;
}
Env
My paraview version is
Client Information:
Version: 5.11.1
VTK Version: 9.2.20220823
Qt Version: 5.15.2
vtkIdType size: 64bits
Embedded Python: On
Python Library Path: /public/software/ParaView-5.11.1/lib/python3.9
Python Library Version: 3.9.13 (main, Mar 30 2023, 16:46:40) [GCC 7.3.1 20180303 (Red Hat 7.3.1-5)]
Python Numpy Support: On
Python Numpy Path: /public/software/ParaView-5.11.1/lib/python3.9/site-packages/numpy
Python Numpy Version: 1.21.1
Python Matplotlib Support: On
Python Matplotlib Path: /public/software/ParaView-5.11.1/lib/python3.9/site-packages/matplotlib
Python Matplotlib Version: 3.2.1
Python Testing: Off
MPI Enabled: On
ParaView Build ID: superbuild 2c64f06366aa320e98d98577af18bad64ad9a630 (!1076)
Disable Registry: Off
Test Directory:
Data Directory:
SMP Backend: TBB
SMP Max Number of Threads: 52
OpenGL Vendor: VMware, Inc.
OpenGL Version: 3.3 (Core Profile) Mesa 18.0.5
OpenGL Renderer: llvmpipe (LLVM 6.0, 256 bits)
Accelerated filters overrides available: No
Connection Information:
Remote Connection: No
The VTK version of my solver links to is 9.1.0
The pyevtk version I used is 1.5.0
File
created_by_VTK.vts (86.9 KB)
created_by_pyevtk.vts (110.0 KB)