A 3D graphics GUI!
(But fewer dimensions ok too!)
The set of programming tools built around Paraview to achieve scientific goals: VTK-m, Fides, ADIOS.
Paraview is the crown jewel of the DOE software ecosystem.
It shows we can meet our unique challenges with performant and usable software.
Easy (on Ubuntu)
$ sudo apt-get install qt5-default libqt5x11extras5-dev \
libqt5help5 qttools5-dev qtxmlpatterns5-dev-tools libqt5svg5-dev
$ git clone --recursive https://gitlab.kitware.com/paraview/paraview.git
$ cd paraview; mkdir build
paraview/build$ cmake ../ -G Ninja
paraview/build$ ninja
paraview/build$ ./bin/paraview
$ git clone https://github.com/NAThompson/paraview_mwes.git
$ cd paraview_mwes
$ git lfs install
$ git lfs pull
$ cd paraview_mwes/example_data
$ ls
ls
1dgraph.vtk euler_spiral.vtk jacobi_theta.vtk scalogram.vtk
daubechies_5_scaling_convergence.csv fermi_surface.vtk lorenz.vtk
- Choose a different colortable
- Color by time, then by curvature
- Use the wireframe representation, and change the linewidth
- Change the background color
- Volume render the
$$E = E(k_x, k_y, k_z)$$ data - Extract a Fermi surface from band21 via isocontouring
- Color the Fermi surface by the gradients
- Extract another Fermi surface from band23, and display it simultaneously with the Fermi surface from band21.
- Label the axes as
$$k_x, k_y, k_z$$ , and color table$$|\nabla E|$$ using$$\LaTeX$$ . - Visualize on the Looking Glass
- Triangulate unstructured points
- Use raytracing to improve image quality
Idiom: Take a source, and apply a filter, create a render.
Yes, we used canned data packaged with the software.
Let's upload some real data.
Easiest thing: A .csv
:
➜ paraview_tutorial$ cat daubechies_5_scaling_convergence.csv | more
r, linear, quadratic_b_spline, cubic_b_spline, quintic_b_spline, cubic_hermite, pchip, makima, quintic_hermite, septic_hermite
2, 0.045726192045372316, 0.024828473821310204, 0.023277651982070491, 0.020662363191588429, 0.003992590683903730, 0.022484271673629985, 0.029204046982060805, 0.010163942298858419, 0.018448258714584442
3, 0.016022890862746997, 0.005755604628244537, 0.005230814883100732, 0.004517849453574918, 0.001041719154716125, 0.012099678303480466, 0.007022457175191787, 0.002578238856418169, 0.004660243256272567
4, 0.004326423048900463, 0.001479662333354892, 0.001403577051129257, 0.001258471400971933, 0.000264701848608107, 0.003008722563990818, 0.002047077677412190, 0.000660995684770599, 0.001198057778804051
5, 0.001309431207276557, 0.000373346431293220, 0.000339254208968076, 0.000291630399409548, 0.000067783984772363, 0.000465714198778056, 0.000491306619344600, 0.000168562863492538, 0.000305161839582041
➜ paraview_tutorial$ cat data/1dgraph.vtk
# vtk DataFile Version 3.0
vtk output
ASCII
DATASET STRUCTURED_GRID
DIMENSIONS 96 1 1 POINTS 96 float
0 0 0
0.0208333 0 0
0.0416667 0 0
...
1.97917 0 0
POINT_DATA 96
SCALARS u double 1
LOOKUP_TABLE default
0.84768
0.787256
...
0.575502
Make a series of .vtk
files, and then label them with an increasing index foo_1.vtk
, foo_3.vtk
, foo_5.vtk
.
Then Paraview knows that this is a dataset that evolves in time.
Let's use this as an example of how to watch the evolution of a simulation in Paraview.
Discretize KdV via Zabusky & Krustal
where
Given a wavelet transform
$$ \mathcal{W}f := \frac{1}{\sqrt{|s|}} \int_{-\infty}^{\infty} f(x) \psi\left( \frac{x-t}{s} \right)^{*} , \mathrm{d}x $$
let's compute a scalogram, a graph of
# vtk DataFile Version 3.0
vtk output
ASCII
DATASET STRUCTURED_GRID
DIMENSIONS 2 2 1 POINTS 4 float
-2 1e-05 0
2 1e-05 0
-2 1 0
2 1 0
POINT_DATA 4
SCALARS scalogram double 1
LOOKUP_TABLE default
6.86438e-06
6.86117e-06
0.000390553
0.000382475
2D image with a field on it can be represented by a colortable, or by a surface.
Let's see how to create a surface.
Let's look at the Jacobi
We might want to look at magnitude, phase, real part, or imaginary part of this function, for various values of
vtkm::cont::DataSet dataSet = ...
// ...
dataSet.AddPointField("r", r.data(), r.size());
dataSet.AddPointField("theta", theta.data(), theta.size());
dataSet.AddPointField("Re(z)", real_part.data(), real_part.size());
dataSet.AddPointField("Im(z)", imag_part.data(), imag_part.size());
Unstructured gets hard in a hurry. Let's solve an easy 1D problem to wrap our heads around it.
Easy to validate: As
Let
This is a dense linear system
double lambda = 2.5;
auto gauss_data = boost::math::quadrature::gauss<double, 256>();
std::vector<double> nodes = ... ;
std::vector<double> weights = ... ;
Eigen::VectorXd f(nodes.size());
for (size_t i = 0; i < nodes.size(); ++i) {
f[i] = rhs(nodes[i]);
}
Eigen::MatrixXd M(nodes.size(), nodes.size());
for (size_t i = 0; i < nodes.size(); ++i) {
for (size_t j = 0; j < nodes.size(); ++j) {
if (i == j) {
M(i,j) = 1;
}
else {
M(i,j) = -lambda*(nodes[i] - nodes[j])*weights[j];
}
}
}
// Solution y at Gauss nodes:
Eigen::VectorXd y = M.fullPivLu().solve(f);
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
for (size_t i = 0; i < nodes.size(); ++i) {
vtkm::Id pid = dsb.AddPoint({nodes[i], 0.0, 0.0});
ids.push_back(pid);
}
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet dataSet;
dataSet = dsb.Create();
dataSet.AddPointField("y", y.data(), y.size());
vtkm::io::VTKDataSetWriter writer("fredholm.vtk");
writer.WriteDataSet(dataSet);
Let's build an Euler spiral:
#include <boost/math/quadrature/tanh_sinh.hpp>
auto integrator = boost::math::quadrature::tanh_sinh<double>();
auto x_coord = [](double s) { return std::cos(s*s); };
auto y_coord = [](double s) { return std::sin(s*s); };
// ...
double x = integrator.integrate(x_coord, 0.0, t);
double y = integrator.integrate(y_coord, 0.0, t);
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
for (size_t i = 0; i < spiral_data.size(); ++i) {
auto point = spiral_data[i];
vtkm::Id pid = dsb.AddPoint({point[0], point[1], 0.0});
ids.push_back(pid);
}
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet dataSet = dsb.Create();
vtkm::io::VTKDataSetWriter writer("euler_spiral.vtk");
writer.WriteDataSet(dataSet);
2D video: The Gray-Scott Initial Value Problem
Subject to
Stability condition
Perfect use case for watching you sim.
Eigen::MatrixXd U0(n,n);
Eigen::MatrixXd U1(n,n);
Eigen::MatrixXd V0(n,n);
Eigen::MatrixXd V1(n,n);
// ...
for (int64_t i = 1; i < n - 1; ++i)
{
for (int64_t j = 1; j < n - 1; ++j)
{
rhsU = Du*(U0(i+1,j) + U0(i,j+1) - 4*U0(i,j) + U0(i-1, j) + U0(i, j-1))/(dx*dx) - U0(i,j)*V0(i,j)*V0(i,j) + F*(1-U0(i,j));
rhsV = Dv*(V0(i+1,j) + V0(i,j+1) - 4*V0(i,j) + V0(i-1, j) + V0(i, j-1))/(dx*dx) + U0(i,j)*V0(i,j)*V0(i,j) - (F+k)*V0(i,j);
U1(i,j) = U0(i,j) + dt*rhsU;
V1(i,j) = V0(i,j) + dt*rhsV;
}
}
vtkm::cont::DataSetBuilderUniform dsb;
vtkm::Id2 dims(n, n);
vtkm::Vec2f_64 origin(0, 0);
double dx = gs_params.x_max/n;
vtkm::Vec2f_64 spacing(dx, dx);
vtkm::cont::DataSet dataSet = dsb.Create(dims, origin, spacing);
dataSet.AddPointField("U", U.data(), U.size());
dataSet.AddPointField("V", V.data(), V.size());
vtkm::io::VTKDataSetWriter writer("gray_scott_" + std::to_string(step_index) + ".vtk");
writer.WriteDataSet(dataSet);
Very simple VTK-m as no connectivity info is needed:
std::random_device rd;
std::normal_distribution<double> dis(0.0, 1.0);
std::array<double, 3> center{0.0, 0.0, 0.0};
vtkm::cont::DataSetBuilderExplicitIterative dsb;
for (size_t i = 0; i < 500; ++i)
{
double x = center[0] + dis(rd);
double y = center[1] + dis(rd);
double z = center[2] + dis(rd);
dsb.AddPoint({x, y, z});
}
vtkm::cont::DataSet dataSet = dsb.Create();
vtkm::io::VTKDataSetWriter writer("scattered.vtk");
writer.WriteDataSet(dataSet);
Start easy with data provided by Paraview.
Further into the VTK-m data model: Rectilinear grids
std::random_device rd;
std::uniform_real_distribution<double> dis(0.0, 1.0/samples);
std::vector<double> x(samples);
std::vector<double> y(samples);
std::vector<double> z(samples);
for (size_t i = 0; i < samples; ++i)
{
x[i] = dis(rd);
y[i] = dis(rd);
z[i] = 0;
}
std::sort(x.begin(), x.end());
std::sort(y.begin(), y.end());
vtkm::cont::DataSetBuilderRectilinear dsb;
vtkm::cont::DataSet ds = dsb.Create(x, y, z);
vtkm::io::VTKDataSetWriter writer("random_rectilinear.vtk");
writer.WriteDataSet(ds);
The Padua points are the union of two Chebyshev grids.
Chebyshev nodes:
$$ \mathrm{Pad}{n} := (C{n+1}^{\mathrm{odd}} \times C_{n+2}^{\mathrm{even}}) \cup (C_{n+1}^{\mathrm{even}} \times C_{n+2}^{\mathrm{odd}}) $$
All points lie on the Lissajous curve
Two Rectilinear
grids; perfect use case for PartitionedDataSet
.
vtkm::cont::DataSetBuilderRectilinear dsb;
std::vector<double> x1((n+1)/2, std::numeric_limits<double>::quiet_NaN());
std::vector<double> y1((n+2)/2, std::numeric_limits<double>::quiet_NaN());
std::vector<double> z1(1, 0);
std::vector<double> x2((n+2)/2, std::numeric_limits<double>::quiet_NaN());
std::vector<double> y2((n+1)/2, std::numeric_limits<double>::quiet_NaN());
std::vector<double> z2(1, 0);
// initialize with Padua points . . .
vtkm::cont::DataSet ds1 = dsb.Create(x1, y1, z1);
vtkm::cont::DataSet ds2 = dsb.Create(x2, y2, z2);
vtkm::cont::PartitionedDataSet pds;
pds.AppendPartitions({ds1, ds2});
for (int i = 0; i < pds.GetNumberOfPartitions(); ++i)
{
auto & p = pds.GetPartition(i);
vtkm::io::VTKDataSetWriter writer("padua_" + std::to_string(i+1) + ".vtk");
writer.WriteDataSet(p);
}
In fluid dynamics, it is common to store pressure as a cell variable, and velocity at a point variable.
Hence we have field associations in VTK-m.
vtkm::cont::Field::Association::POINTS
vtkm::cont::Field::Association::CELL_SET
The .vtk
files are nice, intuitive ASCII.
But they quickly overwhelm the Paraview parser, so we need a binary format to move forward.
vtkm::rendering::Actor actor(dataSet.GetCellSet(),
dataSet.GetCoordinateSystem(),
dataSet.GetField("U"));
vtkm::rendering::Scene scene;
scene.AddActor(actor);
vtkm::rendering::MapperRayTracer mapper;
vtkm::rendering::CanvasRayTracer canvas(1920, 1080);
vtkm::rendering::View2D view(scene, mapper, canvas);
view.Initialize();
view.Paint();
view.SaveAs("gray_scott_u_" + std::to_string(k) + ".pnm");
VTK-m rendering is for sanity checks; it isn't a full-featured renderer like Paraview.
In addition, customizing the graphics in C++ is incredibly painful.