Commit e95f7a60 authored by Petr Strakos's avatar Petr Strakos

adding sources for a simple simulation code, setting user

parent a54de085
cmake_minimum_required(VERSION 3.3)
project(CatalystCxxFullExample)
include_directories(/apps/all/ParaView/5.6.0-intel-2017a-mpi/include/paraview-5.6/)
link_directories(/apps/all/ParaView/5.6.0-intel-2017a-mpi/lib/)
set(USE_CATALYST ON CACHE BOOL "Link the simulator with Catalyst")
if(USE_CATALYST)
find_package(ParaView 4.1 REQUIRED COMPONENTS vtkPVPythonCatalyst)
include("${PARAVIEW_USE_FILE}")
set(Adaptor_SRCS
FEAdaptor.cxx
)
add_library(CxxFullExampleAdaptor ${Adaptor_SRCS})
target_link_libraries(CxxFullExampleAdaptor vtkPVPythonCatalyst vtkParallelMPI)
add_definitions("-DUSE_CATALYST")
if(NOT PARAVIEW_USE_MPI)
message(SEND_ERROR "ParaView must be built with MPI enabled")
endif()
else()
find_package(MPI REQUIRED)
include_directories(${MPI_C_INCLUDE_PATH})
endif()
add_executable(CxxFullExample FEDriver.cxx FEDataStructures.cxx)
if(USE_CATALYST)
target_link_libraries(CxxFullExample LINK_PRIVATE CxxFullExampleAdaptor)
include(vtkModuleMacros)
include(vtkMPI)
vtk_mpi_link(CxxFullExample)
else()
target_link_libraries(CxxFullExample LINK_PRIVATE ${MPI_LIBRARIES})
endif()
option(BUILD_TESTING "Build Testing" OFF)
# Setup testing.
if (BUILD_TESTING)
include(CTest)
add_test(NAME CxxFullExampleTest COMMAND CxxFullExample ${CMAKE_CURRENT_SOURCE_DIR}/SampleScripts/feslicescript.py)
set_tests_properties(CxxFullExampleTest PROPERTIES LABELS "PARAVIEW;CATALYST")
endif()
#include "FEAdaptor.h"
#include "FEDataStructures.h"
#include <iostream>
#include <vtkCPDataDescription.h>
#include <vtkCPInputDataDescription.h>
#include <vtkCPProcessor.h>
#include <vtkCPPythonScriptPipeline.h>
#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkDoubleArray.h>
#include <vtkFloatArray.h>
#include <vtkNew.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkUnstructuredGrid.h>
namespace
{
vtkCPProcessor* Processor = NULL;
vtkUnstructuredGrid* VTKGrid;
void BuildVTKGrid(Grid& grid)
{
// create the points information
vtkNew<vtkDoubleArray> pointArray;
pointArray->SetNumberOfComponents(3);
pointArray->SetArray(
grid.GetPointsArray(), static_cast<vtkIdType>(grid.GetNumberOfPoints() * 3), 1);
vtkNew<vtkPoints> points;
points->SetData(pointArray.GetPointer());
VTKGrid->SetPoints(points.GetPointer());
// create the cells
size_t numCells = grid.GetNumberOfCells();
VTKGrid->Allocate(static_cast<vtkIdType>(numCells * 9));
for (size_t cell = 0; cell < numCells; cell++)
{
unsigned int* cellPoints = grid.GetCellPoints(cell);
vtkIdType tmp[8] = { cellPoints[0], cellPoints[1], cellPoints[2], cellPoints[3], cellPoints[4],
cellPoints[5], cellPoints[6], cellPoints[7] };
VTKGrid->InsertNextCell(VTK_HEXAHEDRON, 8, tmp);
}
}
void UpdateVTKAttributes(Grid& grid, Attributes& attributes, vtkCPInputDataDescription* idd)
{
if (idd->IsFieldNeeded("velocity", vtkDataObject::POINT) == true)
{
if (VTKGrid->GetPointData()->GetNumberOfArrays() == 0)
{
// velocity array
vtkNew<vtkDoubleArray> velocity;
velocity->SetName("velocity");
velocity->SetNumberOfComponents(3);
velocity->SetNumberOfTuples(static_cast<vtkIdType>(grid.GetNumberOfPoints()));
VTKGrid->GetPointData()->AddArray(velocity.GetPointer());
}
vtkDoubleArray* velocity =
vtkDoubleArray::SafeDownCast(VTKGrid->GetPointData()->GetArray("velocity"));
// The velocity array is ordered as vx0,vx1,vx2,..,vy0,vy1,vy2,..,vz0,vz1,vz2,..
// so we need to create a full copy of it with VTK's ordering of
// vx0,vy0,vz0,vx1,vy1,vz1,..
double* velocityData = attributes.GetVelocityArray();
vtkIdType numTuples = velocity->GetNumberOfTuples();
for (vtkIdType i = 0; i < numTuples; i++)
{
double values[3] = { velocityData[i], velocityData[i + numTuples],
velocityData[i + 2 * numTuples] };
velocity->SetTypedTuple(i, values);
}
}
if (idd->IsFieldNeeded("pressure", vtkDataObject::CELL) == true)
{
if (VTKGrid->GetCellData()->GetNumberOfArrays() == 0)
{
// pressure array
vtkNew<vtkFloatArray> pressure;
pressure->SetName("pressure");
pressure->SetNumberOfComponents(1);
VTKGrid->GetCellData()->AddArray(pressure.GetPointer());
}
vtkFloatArray* pressure =
vtkFloatArray::SafeDownCast(VTKGrid->GetCellData()->GetArray("pressure"));
// The pressure array is a scalar array so we can reuse
// memory as long as we ordered the points properly.
float* pressureData = attributes.GetPressureArray();
pressure->SetArray(pressureData, static_cast<vtkIdType>(grid.GetNumberOfCells()), 1);
}
}
void BuildVTKDataStructures(Grid& grid, Attributes& attributes, vtkCPInputDataDescription* idd)
{
if (VTKGrid == NULL)
{
// The grid structure isn't changing so we only build it
// the first time it's needed. If we needed the memory
// we could delete it and rebuild as necessary.
VTKGrid = vtkUnstructuredGrid::New();
BuildVTKGrid(grid);
}
UpdateVTKAttributes(grid, attributes, idd);
}
}
namespace FEAdaptor
{
void Initialize(int numScripts, char* scripts[])
{
if (Processor == NULL)
{
Processor = vtkCPProcessor::New();
Processor->Initialize();
}
else
{
Processor->RemoveAllPipelines();
}
for (int i = 0; i < numScripts; i++)
{
vtkNew<vtkCPPythonScriptPipeline> pipeline;
pipeline->Initialize(scripts[i]);
Processor->AddPipeline(pipeline.GetPointer());
}
}
void Finalize()
{
if (Processor)
{
Processor->Delete();
Processor = NULL;
}
if (VTKGrid)
{
VTKGrid->Delete();
VTKGrid = NULL;
}
}
void CoProcess(
Grid& grid, Attributes& attributes, double time, unsigned int timeStep, bool lastTimeStep)
{
vtkNew<vtkCPDataDescription> dataDescription;
dataDescription->AddInput("input");
dataDescription->SetTimeData(time, timeStep);
if (lastTimeStep == true)
{
// assume that we want to all the pipelines to execute if it
// is the last time step.
dataDescription->ForceOutputOn();
}
if (Processor->RequestDataDescription(dataDescription.GetPointer()) != 0)
{
vtkCPInputDataDescription* idd = dataDescription->GetInputDescriptionByName("input");
BuildVTKDataStructures(grid, attributes, idd);
idd->SetGrid(VTKGrid);
Processor->CoProcess(dataDescription.GetPointer());
}
}
} // end of Catalyst namespace
#ifndef FEADAPTOR_HEADER
#define FEADAPTOR_HEADER
class Attributes;
class Grid;
namespace FEAdaptor
{
void Initialize(int numScripts, char* scripts[]);
void Finalize();
void CoProcess(
Grid& grid, Attributes& attributes, double time, unsigned int timeStep, bool lastTimeStep);
}
#endif
#include "FEDataStructures.h"
#include <iostream>
#include <mpi.h>
Grid::Grid()
{
}
void Grid::Initialize(const unsigned int numPoints[3], const double spacing[3])
{
if (numPoints[0] == 0 || numPoints[1] == 0 || numPoints[2] == 0)
{
std::cerr << "Must have a non-zero amount of points in each direction.\n";
}
// in parallel, we do a simple partitioning in the x-direction.
int mpiSize = 1;
int mpiRank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
unsigned int startXPoint = mpiRank * numPoints[0] / mpiSize;
unsigned int endXPoint = (mpiRank + 1) * numPoints[0] / mpiSize;
if (mpiSize != mpiRank + 1)
{
endXPoint++;
}
// create the points -- slowest in the x and fastest in the z directions
double coord[3] = { 0, 0, 0 };
for (unsigned int i = startXPoint; i < endXPoint; i++)
{
coord[0] = i * spacing[0];
for (unsigned int j = 0; j < numPoints[1]; j++)
{
coord[1] = j * spacing[1];
for (unsigned int k = 0; k < numPoints[2]; k++)
{
coord[2] = k * spacing[2];
// add the coordinate to the end of the vector
std::copy(coord, coord + 3, std::back_inserter(this->Points));
}
}
}
// create the hex cells
unsigned int cellPoints[8];
unsigned int numXPoints = endXPoint - startXPoint;
for (unsigned int i = 0; i < numXPoints - 1; i++)
{
for (unsigned int j = 0; j < numPoints[1] - 1; j++)
{
for (unsigned int k = 0; k < numPoints[2] - 1; k++)
{
cellPoints[0] = i * numPoints[1] * numPoints[2] + j * numPoints[2] + k;
cellPoints[1] = (i + 1) * numPoints[1] * numPoints[2] + j * numPoints[2] + k;
cellPoints[2] = (i + 1) * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k;
cellPoints[3] = i * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k;
cellPoints[4] = i * numPoints[1] * numPoints[2] + j * numPoints[2] + k + 1;
cellPoints[5] = (i + 1) * numPoints[1] * numPoints[2] + j * numPoints[2] + k + 1;
cellPoints[6] = (i + 1) * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k + 1;
cellPoints[7] = i * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k + 1;
std::copy(cellPoints, cellPoints + 8, std::back_inserter(this->Cells));
}
}
}
}
size_t Grid::GetNumberOfPoints()
{
return this->Points.size() / 3;
}
size_t Grid::GetNumberOfCells()
{
return this->Cells.size() / 8;
}
double* Grid::GetPointsArray()
{
if (this->Points.empty())
{
return NULL;
}
return &(this->Points[0]);
}
double* Grid::GetPoint(size_t pointId)
{
if (pointId >= this->Points.size())
{
return NULL;
}
return &(this->Points[pointId * 3]);
}
unsigned int* Grid::GetCellPoints(size_t cellId)
{
if (cellId >= this->Cells.size())
{
return NULL;
}
return &(this->Cells[cellId * 8]);
}
Attributes::Attributes()
{
this->GridPtr = NULL;
}
void Attributes::Initialize(Grid* grid)
{
this->GridPtr = grid;
}
void Attributes::UpdateFields(double time)
{
size_t numPoints = this->GridPtr->GetNumberOfPoints();
this->Velocity.resize(numPoints * 3);
// provide different update setting for different parallel process
int mpiSize = 1;
int mpiRank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
double setting = 1.0 + (double) mpiRank / (double) mpiSize;
for (size_t pt = 0; pt < numPoints; pt++)
{
double* coord = this->GridPtr->GetPoint(pt);
this->Velocity[pt] = coord[1] * time * setting;
}
std::fill(this->Velocity.begin() + numPoints, this->Velocity.end(), 0.0);
size_t numCells = this->GridPtr->GetNumberOfCells();
this->Pressure.resize(numCells);
std::fill(this->Pressure.begin(), this->Pressure.end(), setting);
}
double* Attributes::GetVelocityArray()
{
if (this->Velocity.empty())
{
return NULL;
}
return &this->Velocity[0];
}
float* Attributes::GetPressureArray()
{
if (this->Pressure.empty())
{
return NULL;
}
return &this->Pressure[0];
}
#ifndef FEDATASTRUCTURES_HEADER
#define FEDATASTRUCTURES_HEADER
#include <cstddef>
#include <vector>
class Grid
{
public:
Grid();
void Initialize(const unsigned int numPoints[3], const double spacing[3]);
size_t GetNumberOfPoints();
size_t GetNumberOfCells();
double* GetPointsArray();
double* GetPoint(size_t pointId);
unsigned int* GetCellPoints(size_t cellId);
private:
std::vector<double> Points;
std::vector<unsigned int> Cells;
};
class Attributes
{
// A class for generating and storing point and cell fields.
// Velocity is stored at the points and pressure is stored
// for the cells. The current velocity profile is for a
// shearing flow with U(y,t) = y*t, V = 0 and W = 0.
// Pressure is constant through the domain.
public:
Attributes();
void Initialize(Grid* grid);
void UpdateFields(double time);
double* GetVelocityArray();
float* GetPressureArray();
private:
std::vector<double> Velocity;
std::vector<float> Pressure;
Grid* GridPtr;
};
#endif
#include "FEDataStructures.h"
#include <mpi.h>
#include <stdio.h>
#include <unistd.h>
#include <iostream>
#include <stdlib.h>
#ifdef USE_CATALYST
#include "FEAdaptor.h"
#endif
// Example of a C++ adaptor for a simulation code
int main(int argc, char** argv)
{
// Check the input arguments for area size
if (argc < 4) {
printf("Not all arguments for grid definition supplied\n");
return 0;
}
unsigned int pointsX = abs(std::stoi(argv[1]));
unsigned int pointsY = abs(std::stoi(argv[2]));
unsigned int pointsZ = abs(std::stoi(argv[3]));
//MPI_Init(&argc, &argv);
MPI_Init(NULL, NULL);
Grid grid;
unsigned int numPoints[3] = { pointsX, pointsY, pointsZ };
double spacing[3] = { 1, 1.1, 1.3 };
grid.Initialize(numPoints, spacing);
Attributes attributes;
attributes.Initialize(&grid);
#ifdef USE_CATALYST
// The first argument is the program name
FEAdaptor::Initialize(argc - 4, &argv[4]);
#endif
unsigned int numberOfTimeSteps = 1000;
for (unsigned int timeStep = 0; timeStep < numberOfTimeSteps; timeStep++)
{
// use a time step length of 0.1
double time = timeStep * 0.1;
attributes.UpdateFields(time);
#ifdef USE_CATALYST
FEAdaptor::CoProcess(grid, attributes, time, timeStep, timeStep == numberOfTimeSteps - 1);
#endif
// Get the name of the processor
char processor_name[MPI_MAX_PROCESSOR_NAME];
int name_len;
MPI_Get_processor_name(processor_name, &name_len);
printf("This is processor %s, time step: %0.3f\n", processor_name, time);
usleep(500000);
}
#ifdef USE_CATALYST
FEAdaptor::Finalize();
#endif
MPI_Finalize();
return 0;
}
from paraview.simple import *
from paraview import coprocessing
#--------------------------------------------------------------
# Code generated from cpstate.py to create the CoProcessor.
# ----------------------- CoProcessor definition -----------------------
def CreateCoProcessor():
def _CreatePipeline(coprocessor, datadescription):
class Pipeline:
filename_3_pvtu = coprocessor.CreateProducer( datadescription, "input" )
Slice1 = Slice( guiName="Slice1", Crinkleslice=0, SliceOffsetValues=[0.0], Triangulatetheslice=1, SliceType="Plane" )
Slice1.SliceType.Offset = 0.0
Slice1.SliceType.Origin = [34.5, 32.45, 27.95]
Slice1.SliceType.Normal = [1.0, 0.0, 0.0]
# create a new 'Parallel PolyData Writer'
parallelPolyDataWriter1 = servermanager.writers.XMLPPolyDataWriter(Input=Slice1)
# register the writer with coprocessor
# and provide it with information such as the filename to use,
# how frequently to write the data, etc.
coprocessor.RegisterWriter(parallelPolyDataWriter1, filename='slice_%t.pvtp', freq=10)
# create a new 'Parallel UnstructuredGrid Writer'
unstructuredGridWriter1 = servermanager.writers.XMLPUnstructuredGridWriter(Input=filename_3_pvtu)
# register the writer with coprocessor
# and provide it with information such as the filename to use,
# how frequently to write the data, etc.
coprocessor.RegisterWriter(unstructuredGridWriter1, filename='fullgrid_%t.pvtu', freq=100)
return Pipeline()
class CoProcessor(coprocessing.CoProcessor):
def CreatePipeline(self, datadescription):
self.Pipeline = _CreatePipeline(self, datadescription)
coprocessor = CoProcessor()
freqs = {'input': [10, 100]}
coprocessor.SetUpdateFrequencies(freqs)
return coprocessor
#--------------------------------------------------------------
# Global variables that will hold the pipeline for each timestep
# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
# It will be automatically setup when coprocessor.UpdateProducers() is called the
# first time.
coprocessor = CreateCoProcessor()
#--------------------------------------------------------------
# Enable Live-Visualizaton with ParaView
coprocessor.EnableLiveVisualization(False)
# ---------------------- Data Selection method ----------------------
def RequestDataDescription(datadescription):
"Callback to populate the request for current timestep"
global coprocessor
if datadescription.GetForceOutput() == True:
# We are just going to request all fields and meshes from the simulation
# code/adaptor.
for i in range(datadescription.GetNumberOfInputDescriptions()):
datadescription.GetInputDescription(i).AllFieldsOn()
datadescription.GetInputDescription(i).GenerateMeshOn()
return
# setup requests for all inputs based on the requirements of the
# pipeline.
coprocessor.LoadRequestedData(datadescription)
# ------------------------ Processing method ------------------------
def DoCoProcessing(datadescription):
"Callback to do co-processing for current timestep"
global coprocessor
# Update the coprocessor by providing it the newly generated simulation data.
# If the pipeline hasn't been setup yet, this will setup the pipeline.
coprocessor.UpdateProducers(datadescription)
# Write output data, if appropriate.
coprocessor.WriteData(datadescription);
# Write image capture (Last arg: rescale lookup table), if appropriate.
coprocessor.WriteImages(datadescription, rescale_lookuptable=False)
# Live Visualization, if enabled.
coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
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