diff --git a/docs.it4i/software/viz/insitu/CMakeLists.txt b/docs.it4i/software/viz/insitu/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..c1f7b97dfd678898b20c0000b1e94edb0d2cd7a8 --- /dev/null +++ b/docs.it4i/software/viz/insitu/CMakeLists.txt @@ -0,0 +1,41 @@ +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() diff --git a/docs.it4i/software/viz/insitu/FEAdaptor.cxx b/docs.it4i/software/viz/insitu/FEAdaptor.cxx new file mode 100644 index 0000000000000000000000000000000000000000..6409f18ff87dadd1ea191ed1a3091f387582334d --- /dev/null +++ b/docs.it4i/software/viz/insitu/FEAdaptor.cxx @@ -0,0 +1,162 @@ +#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 diff --git a/docs.it4i/software/viz/insitu/FEAdaptor.h b/docs.it4i/software/viz/insitu/FEAdaptor.h new file mode 100644 index 0000000000000000000000000000000000000000..9bc2770484a55b3f117a3bd04ce1107706e8cd43 --- /dev/null +++ b/docs.it4i/software/viz/insitu/FEAdaptor.h @@ -0,0 +1,17 @@ +#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 diff --git a/docs.it4i/software/viz/insitu/FEDataStructures.cxx b/docs.it4i/software/viz/insitu/FEDataStructures.cxx new file mode 100644 index 0000000000000000000000000000000000000000..858e450aa94dea6ccc598e0387596e2c7b96d15d --- /dev/null +++ b/docs.it4i/software/viz/insitu/FEDataStructures.cxx @@ -0,0 +1,155 @@ +#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]; +} diff --git a/docs.it4i/software/viz/insitu/FEDataStructures.h b/docs.it4i/software/viz/insitu/FEDataStructures.h new file mode 100644 index 0000000000000000000000000000000000000000..b9d39db966fb14c1a3e86c9aa2cd79289ae4f843 --- /dev/null +++ b/docs.it4i/software/viz/insitu/FEDataStructures.h @@ -0,0 +1,42 @@ +#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 diff --git a/docs.it4i/software/viz/insitu/FEDriver.cxx b/docs.it4i/software/viz/insitu/FEDriver.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9c7d9ec4e21770daddae374959381a3b1854b375 --- /dev/null +++ b/docs.it4i/software/viz/insitu/FEDriver.cxx @@ -0,0 +1,65 @@ +#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; +} diff --git a/docs.it4i/software/viz/insitu/feslicescript.py b/docs.it4i/software/viz/insitu/feslicescript.py new file mode 100644 index 0000000000000000000000000000000000000000..21e9d7522339216c1cdd12cb24bef232d7e9f1d5 --- /dev/null +++ b/docs.it4i/software/viz/insitu/feslicescript.py @@ -0,0 +1,93 @@ +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) diff --git a/docs.it4i/software/viz/insitu/insitu.tar.gz b/docs.it4i/software/viz/insitu/insitu.tar.gz new file mode 100644 index 0000000000000000000000000000000000000000..d207465c54760da5dc15d3fae7c05591f7bf5547 Binary files /dev/null and b/docs.it4i/software/viz/insitu/insitu.tar.gz differ