Coprocessing example

This example is used to demonstrate how the co-processing library can be used with a C++ based simulation code. In the ParaView/CoProcessing/Adaptors/FortranAdaptors directory there is code useful for integrating C or Fortran based simulation codes withthe co-processing library. Note that this example requires MPI to be available on your system. The executable takes in a python coprocessing script and a number of time steps to be run for. Note to remember to set your system environment properly. See [[1]] for details.

CoProcessingExample.cxx

<source lang="cpp">

1. include "vtkCPDataDescription.h"
2. include "vtkCPInputDataDescription.h"
3. include "vtkCPProcessor.h"
4. include "vtkCPPythonScriptPipeline.h"
5. include "vtkElevationFilter.h"
6. include "vtkMultiProcessController.h"
7. include "vtkPolyData.h"
8. include "vtkSmartPointer.h"
9. include "vtkSphereSource.h"
10. include "vtkXMLUnstructuredGridReader.h"

1. include <mpi.h>
2. include <string>

class DataGenerator { public:

 DataGenerator()
   {
   this->Sphere = vtkSmartPointer<vtkSphereSource>::New();
   this->Sphere->SetThetaResolution(30);
   this->Sphere->SetPhiResolution(30);
   int procId;
   MPI_Comm_rank(MPI_COMM_WORLD, &procId);
   this->Sphere->SetCenter(procId*4.0, 0, 0);
   this->Elevation = vtkSmartPointer<vtkElevationFilter>::New();
   this->Elevation->SetInputConnection(this->Sphere->GetOutputPort());
   this->Index = 0;
   }

 vtkSmartPointer<vtkPolyData> GetNext()
   {
   double radius = fabs(sin(0.1 * this->Index));
   this->Index++;
   this->Sphere->SetRadius(1.0 + radius);
   this->Elevation->Update();
   vtkSmartPointer<vtkPolyData> ret = vtkSmartPointer<vtkPolyData>::New();
   ret->DeepCopy(this->Elevation->GetOutput());
   return ret;
   }

protected:

 int Index;
 vtkSmartPointer<vtkSphereSource> Sphere;
 vtkSmartPointer<vtkElevationFilter> Elevation;

};

int main(int argc, char* argv[]) {

 if (argc < 3)
   {
   printf("Usage: %s <cp python file> <number of steps>\n", argv[0]);
   return 1;
   }
 // we assume that this is done in parallel
 MPI_Init(&argc, &argv);

 std::string cpPythonFile = argv[1];
 int nSteps = atoi(argv[2]);

 vtkCPProcessor* processor = vtkCPProcessor::New();
 processor->Initialize();
 vtkCPPythonScriptPipeline* pipeline = vtkCPPythonScriptPipeline::New();

 // read the coprocessing python file
 if(pipeline->Initialize(cpPythonFile.c_str()) == 0)
   {
   cout << "Problem reading the python script.\n";
   return 1;
   }

 processor->AddPipeline(pipeline);
 pipeline->Delete();

 if (nSteps == 0)
   {
   return 0;
   }

 // create a data source, typically this will come from the adaptor
 // but here we use generator to create it ourselves
 DataGenerator generator;

 // do coprocessing
 double tStart = 0.0;
 double tEnd = 1.0;
 double stepSize = (tEnd - tStart)/nSteps;

 vtkCPDataDescription* dataDesc = vtkCPDataDescription::New();
 dataDesc->AddInput("input");

 for (int i = 0; i < nSteps; ++i)
   {
   double currentTime = tStart + stepSize*i;
   // set the current time and time step
   dataDesc->SetTimeData(currentTime, i);

   // the last time step we want to force all the pipelines
   // to execute
   if(i == nSteps-1)
     {
     dataDesc->ForceOutputOn();
     }

   // check if the script says we should do coprocessing now
   if(processor->RequestDataDescription(dataDesc) != 0)
     {
     // we are going to do coprocessing so use generator to
     // create our grid at this timestep and provide it to
     // the coprocessing library
     vtkSmartPointer<vtkDataObject> dataObject =
       generator.GetNext();

     dataDesc->GetInputDescriptionByName("input")->SetGrid(dataObject);
     processor->CoProcess(dataDesc);
     }
   }

 dataDesc->Delete();
 processor->Finalize();
 processor->Delete();

 MPI_Finalize();

 return 0;

} </source>

CMakeLists.txt

<source lang="cmake"> cmake_minimum_required(VERSION 2.6)

PROJECT(CoProcessingExample)

FIND_PACKAGE(ParaView REQUIRED) INCLUDE(${PARAVIEW_USE_FILE})

INCLUDE_DIRECTORIES(${VTK_MPI_INCLUDE_DIR}) IF(NOT VTK_MPI_INCLUDE_DIR)

 MESSAGE(SEND_ERROR "Must build ParaView with MPI.")

ENDIF()

ADD_EXECUTABLE(CoProcessingExample CoProcessingExample.cxx) TARGET_LINK_LIBRARIES(CoProcessingExample vtkCoProcessor) </source>

Python Scripts

The first python script below is used to just output the actual results of the example. This would correspond to a simulation run with a coarse grid in order to set up coprocessing runs for larger grids where outputting the entire simulation results can be computationally prohibitive.

<source lang="python"> try: paraview.simple except: from paraview.simple import *

def RequestDataDescription(datadescription):

   "Callback to populate the request for current timestep"
   if datadescription.GetForceOutput() == True:
       for i in range(datadescription.GetNumberOfInputDescriptions()):
           datadescription.GetInputDescription(i).AllFieldsOn()
           datadescription.GetInputDescription(i).GenerateMeshOn()
       return

   timestep = datadescription.GetTimeStep()
   input_name = 'input'
   if (timestep % 1 == 0) :
       datadescription.GetInputDescriptionByName(input_name).AllFieldsOn()
       datadescription.GetInputDescriptionByName(input_name).GenerateMeshOn()
   else:
       datadescription.GetInputDescriptionByName(input_name).AllFieldsOff()
       datadescription.GetInputDescriptionByName(input_name).GenerateMeshOff()

def DoCoProcessing(datadescription):

   "Callback to do co-processing for current timestep"
   cp_writers = []
   timestep = datadescription.GetTimeStep()

   Sphere1 = CreateProducer( datadescription, "input" )
   ParallelPolyDataWriter1 = CreateWriter( XMLPPolyDataWriter, "input_grid_%t.pvtp", 1, cp_writers )

   for writer in cp_writers:
       if timestep % writer.cpFrequency == 0 or datadescription.GetForceOutput() == True:
           writer.FileName = writer.cpFileName.replace("%t", str(timestep))
           writer.UpdatePipeline()

   # explicitly delete the proxies -- we do it this way to avoid problems with prototypes
   tobedeleted = GetNextProxyToDelete()
   while tobedeleted != None:
       Delete(tobedeleted)
       tobedeleted = GetNextProxyToDelete()

def GetNextProxyToDelete():

   iter = servermanager.vtkSMProxyIterator()
   iter.Begin()
   while not iter.IsAtEnd():
     if iter.GetGroup().find("prototypes") != -1:
        iter.Next()
        continue
     proxy = servermanager._getPyProxy(iter.GetProxy())
     proxygroup = iter.GetGroup()
     iter.Next()
     if proxygroup != 'timekeeper' and proxy != None and proxygroup.find("pq_helper_proxies") == -1 :
         return proxy

   return None

def CreateProducer(datadescription, gridname):

 "Creates a producer proxy for the grid"
 if not datadescription.GetInputDescriptionByName(gridname):
   raise RuntimeError, "Simulation input name '%s' does not exist" % gridname
 grid = datadescription.GetInputDescriptionByName(gridname).GetGrid()
 producer = TrivialProducer()
 producer.GetClientSideObject().SetOutput(grid)
 producer.UpdatePipeline()
 return producer

def CreateWriter(proxy_ctor, filename, freq, cp_writers):

   writer = proxy_ctor()
   writer.FileName = filename
   writer.add_attribute("cpFrequency", freq)
   writer.add_attribute("cpFileName", filename)
   cp_writers.append(writer)
   return writer

</source>

---

This second script is still rather simple and only performs a cut on the input from the simulation code. It demonstrates though how desired results can be obtained while performing coprocessing at specified time steps.

<source lang="python"> try: paraview.simple except: from paraview.simple import *

def RequestDataDescription(datadescription):

   "Callback to populate the request for current timestep"
   if datadescription.GetForceOutput() == True:
       for i in range(datadescription.GetNumberOfInputDescriptions()):
           datadescription.GetInputDescription(i).AllFieldsOn()
           datadescription.GetInputDescription(i).GenerateMeshOn()
       return

   timestep = datadescription.GetTimeStep()
   input_name = 'input'
   if (timestep % 5 == 0) :
       datadescription.GetInputDescriptionByName(input_name).AllFieldsOn()
       datadescription.GetInputDescriptionByName(input_name).GenerateMeshOn()
   else:
       datadescription.GetInputDescriptionByName(input_name).AllFieldsOff()
       datadescription.GetInputDescriptionByName(input_name).GenerateMeshOff()

def DoCoProcessing(datadescription):

   "Callback to do co-processing for current timestep"
   cp_writers = []
   timestep = datadescription.GetTimeStep()

   filename_0_pvtp = CreateProducer( datadescription, "input" )

   Clip2 = Clip( guiName="Clip2", InsideOut=0, UseValueAsOffset=0, Scalars=['POINTS', 'Elevation'], Value=0.0, ClipType="Plane" )
   Clip2.ClipType.Normal = [0.0, 1.0, 0.0]
   Clip2.ClipType.Origin = [1.9999999105930328, 0.0, 0.0]
   Clip2.ClipType.Offset = 0.0

   ParallelUnstructuredGridWriter2 = CreateWriter( XMLPUnstructuredGridWriter, "Cut_%t.pvtu", 5, cp_writers )
   for writer in cp_writers:
       if timestep % writer.cpFrequency == 0 or datadescription.GetForceOutput() == True:
           writer.FileName = writer.cpFileName.replace("%t", str(timestep))
           writer.UpdatePipeline()

   # explicitly delete the proxies -- we do it this way to avoid problems with prototypes
   tobedeleted = GetNextProxyToDelete()
   while tobedeleted != None:
       Delete(tobedeleted)
       tobedeleted = GetNextProxyToDelete()

def GetNextProxyToDelete():

   iter = servermanager.vtkSMProxyIterator()
   iter.Begin()
   while not iter.IsAtEnd():
     if iter.GetGroup().find("prototypes") != -1:
        iter.Next()
        continue
     proxy = servermanager._getPyProxy(iter.GetProxy())
     proxygroup = iter.GetGroup()
     iter.Next()
     if proxygroup != 'timekeeper' and proxy != None and proxygroup.find("pq_helper_proxies") == -1 :
         return proxy

   return None

def CreateProducer(datadescription, gridname):

 "Creates a producer proxy for the grid"
 if not datadescription.GetInputDescriptionByName(gridname):
   raise RuntimeError, "Simulation input name '%s' does not exist" % gridname
 grid = datadescription.GetInputDescriptionByName(gridname).GetGrid()
 producer = TrivialProducer()
 producer.GetClientSideObject().SetOutput(grid)
 producer.UpdatePipeline()
 return producer

def CreateWriter(proxy_ctor, filename, freq, cp_writers):

   writer = proxy_ctor()
   writer.FileName = filename
   writer.add_attribute("cpFrequency", freq)
   writer.add_attribute("cpFileName", filename)
   cp_writers.append(writer)
   return writer

</source>