Commit 337a0caf authored by Robert Maynard's avatar Robert Maynard
Browse files

Adding the VTKm Accelerators module.

parent 1e3e2f69
##=============================================================================
##
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
##
## This software is distributed WITHOUT ANY WARRANTY; without even
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the above copyright notice for more information.
##
## Copyright 2012 Sandia Corporation.
## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
## the U.S. Government retains certain rights in this software.
##
##=============================================================================
#ensure we link against our dependencies
include(module.cmake)
find_package(VTKm REQUIRED
OPTIONAL_COMPONENTS Serial CUDA TBB
)
set(lib_srcs
vtkmlib/PolyDataConverter.cxx
vtkmlib/UnstructuredGridConverter.cxx
vtkmlib/ArrayConverters.cxx
vtkmlib/CellSetConverters.cxx
vtkmlib/DataSetConverters.cxx
vtkmlib/Storage.cxx
)
#needed to properly setup language wrappers
set(headers
vtkmContour.h
vtkmThreshold.h
vtkmLevelOfDetail.h
vtkmCellAverage.h
vtkmGradient.h
)
#implementation of the algorithms for cpu accelerators
set(cpu_accelerator_srcs
vtkmContour.cxx
vtkmThreshold.cxx
vtkmLevelOfDetail.cxx
vtkmCellAverage.cxx
vtkmCellSetExplicit.cxx
vtkmCellSetSingleType.cxx
vtkmConnectivityExec.cxx
vtkmGradient.cxx
vtkmlib/Portals.cxx
)
#implementation of the algorithms for gpu accelerators
set(cuda_accelerator_srcs
vtkmContour.cu
vtkmThreshold.cu
vtkmLevelOfDetail.cu
vtkmCellAverage.cu
vtkmCellSetExplicit.cu
vtkmCellSetSingleType.cu
vtkmConnectivityExec.cu
vtkmGradient.cu
vtkmlib/Portals.cu
)
set(VTKM_FILTER_INCLUDE_AOS ${VTK_DISPATCH_AOS_ARRAYS})
set(VTKM_FILTER_INCLUDE_SOA ${VTK_DISPATCH_SOA_ARRAYS})
set(VTKM_FILTERS_ARE_BUILT_STATIC false)
if(BUILD_SHARED_LIBS)
set(VTKM_FILTERS_ARE_BUILT_STATIC true)
endif()
configure_file("${CMAKE_CURRENT_SOURCE_DIR}/vtkmConfig.h.in"
"${CMAKE_CURRENT_BINARY_DIR}/vtkmConfig.h" @ONLY)
#mark all the helper classes as being excluded from wrappers
set_source_files_properties(
vtkmlib/PolyDataConverter
vtkmlib/UnstructuredGridConverter
vtkmlib/ArrayConverters
vtkmlib/CellSetConverters
vtkmlib/DataSetConverters
vtkmlib/Storage
vtkmlib/Portals
vtkmCellSetExplicit
vtkmCellSetSingleType
vtkmConnectivityExec
PROPERTIES
WRAP_EXCLUDE 1
WRAP_EXCLUDE_PYTHON 1
)
set(${vtk-module}_HDRS
vtkmTags.h
vtkmFilterPolicy.h
${CMAKE_CURRENT_BINARY_DIR}/vtkmConfig.h
)
#we are building with CUDA support
if(VTKm_CUDA_FOUND)
#need to find cudadevrt
find_library(CUDA_cudadevrt_LIBRARY cudadevrt
PATHS ${CUDA_TOOLKIT_TARGET_DIR}
PATH_SUFFIXES "x64" "lib64" "libx64"
)
########
## cache and clear the CUDA_NVCC_FLAGS so that they aren't passed to
## the linker. FINDCUDA has some problems properly unquoting CUDA_NVCC_FLAGS
## when "generate-code arch..." is used, so we instead patch the options
##
########
set(compile_options)
foreach(f ${CUDA_NVCC_FLAGS})
if(f MATCHES "generate-code ")
string(REPLACE "generate-code " "generate-code=" f "${f}")
endif()
list(APPEND compile_options ${f})
endforeach()
if(BUILD_SHARED_LIBS AND NOT WIN32)
list(APPEND compile_options -Xcompiler=${CMAKE_CXX_COMPILE_OPTIONS_VISIBILITY}hidden)
list(APPEND compile_options -Xcompiler=-fPIC)
#nvcc doesn't like the macros in VTK and generates hundreds of warnings
#that are false positives
list(APPEND compile_options --disable-warnings)
endif()
set(seperable_state ${CUDA_SEPARABLE_COMPILATION})
set(cache_flag_state ${CUDA_NVCC_FLAGS})
set(CUDA_NVCC_FLAGS "")
set(CUDA_SEPARABLE_COMPILATION ON)
#Some versions of VTK-m overload the CUDA_LIBRARIES to contain private
if(PRIVATE IN_LIST CUDA_LIBRARIES)
set(cache_cuda_libs ${CUDA_LIBRARIES})
set(cache_devrt_libs ${CUDA_cudadevrt_LIBRARY})
set(CUDA_LIBRARIES ${CUDA_LIBRARIES} vtkFiltersGeneral)
set(CUDA_cudadevrt_LIBRARY PRIVATE ${CUDA_cudadevrt_LIBRARY})
endif()
# CUDA doesn't obey usage requirements so we have to use
# CUDA_INCLUDE_DIRECTORIES, but do get the proper list of
# include dirs I need to query the module system, which
# doesnt exist currently, so we manually call vtk_module_impl
vtk_module_impl()
cuda_include_directories(${CMAKE_CURRENT_BINARY_DIR}
${CMAKE_CURRENT_SOURCE_DIR}
${VTKm_INCLUDE_DIRS}
${vtkAcceleratorsVTKm_DEPENDS_INCLUDE_DIRS})
cuda_add_library(vtkAcceleratorsVTKmCuda STATIC
${cuda_accelerator_srcs}
OPTIONS "${compile_options}"
)
set_target_properties(vtkAcceleratorsVTKmCuda
PROPERTIES POSITION_INDEPENDENT_CODE True)
vtk_module_library(vtkAcceleratorsVTKm
${headers}
${lib_srcs}
)
target_link_libraries(vtkAcceleratorsVTKm
PRIVATE vtkAcceleratorsVTKmCuda ${cache_devrt_libs})
set(CUDA_SEPARABLE_COMPILATION ${seperable_state})
set(CUDA_NVCC_FLAGS_CACHE ${cache_flag_state})
if(cache_cuda_libs)
set(CUDA_LIBRARIES ${cache_cuda_libs})
set(CUDA_cudadevrt_LIBRARY ${CUDA_cudadevrt_LIBRARY})
endif()
else()
vtk_module_library(vtkAcceleratorsVTKm
${headers}
${lib_srcs}
${cpu_accelerator_srcs}
)
endif()
#We need to system up VTK-m as a system include dir so that modules
#such as wrapping that depend on vtkAcceleratorsVTKm properly find
#the headers
target_include_directories(vtkAcceleratorsVTKm PUBLIC ${VTKm_INCLUDE_DIRS})
vtk_module_link_libraries(vtkAcceleratorsVTKm LINK_PRIVATE ${VTKm_LIBRARIES})
target_compile_options(vtkAcceleratorsVTKm PRIVATE ${VTKm_COMPILE_OPTIONS})
#install the required headers to make your own vtkm-vtk filter
if(NOT VTK_INSTALL_NO_DEVELOPMENT)
install(DIRECTORY
${CMAKE_CURRENT_SOURCE_DIR}/vtkmlib
DESTINATION ${VTK_INSTALL_INCLUDE_DIR}
COMPONENT Development)
endif()
find_package(VTKm REQUIRED COMPONENTS Base)
include_directories(${VTKm_INCLUDE_DIRS})
vtk_add_test_cxx(${vtk-module}CxxTests tests
# TestVTKMGradientAndVorticity.cxx,NO_VALID
TestVTKMLevelOfDetail.cxx
TestVTKMMarchingCubes.cxx
TestVTKMThreshold.cxx
TestVTKMThreshold2.cxx
)
vtk_test_cxx_executable(${vtk-module}CxxTests tests
RENDERING_FACTORY
)
/*=========================================================================
Program: Visualization Toolkit
Module: TestGradientAndVorticity.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
/*----------------------------------------------------------------------------
Copyright (c) Sandia Corporation
See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/
#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkDoubleArray.h"
#include "vtkmGradient.h"
#include "vtkPointData.h"
#include "vtkSmartPointer.h"
#include "vtkStdString.h"
#include "vtkStructuredGrid.h"
#include "vtkStructuredGridReader.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridReader.h"
#include <vector>
#define VTK_CREATE(type, var) \
vtkSmartPointer<type> var = vtkSmartPointer<type>::New()
namespace
{
double Tolerance = 0.00001;
bool ArePointsWithinTolerance(double v1, double v2)
{
if(v1 == v2 || fabs(v1)+fabs(v2) < Tolerance)
{
return true;
}
if(v1 == 0.0)
{
if(fabs(v2) < Tolerance)
{
return true;
}
std::cout << fabs(v2) << " (fabs(v2)) should be less than "
<< Tolerance << std::endl;
return false;
}
if(fabs(v1/v2) < Tolerance)
{
return true;
}
std::cout << fabs(v1/v2) << " (fabs(v1/v2)) should be less than "
<< Tolerance << std::endl;
return false;
}
//-----------------------------------------------------------------------------
void CreateCellData(vtkDataSet* grid, int numberOfComponents, int offset,
const char* arrayName)
{
vtkIdType numberOfCells = grid->GetNumberOfCells();
VTK_CREATE(vtkDoubleArray, array);
array->SetNumberOfComponents(numberOfComponents);
array->SetNumberOfTuples(numberOfCells);
std::vector<double> tupleValues(numberOfComponents);
double point[3], parametricCenter[3], weights[100];
for(vtkIdType i=0;i<numberOfCells;i++)
{
vtkCell* cell = grid->GetCell(i);
cell->GetParametricCenter(parametricCenter);
int subId = 0;
cell->EvaluateLocation(subId, parametricCenter, point, weights);
for(int j=0;j<numberOfComponents;j++)
{// +offset makes the curl/vorticity nonzero
tupleValues[j] = point[(j+offset)%3];
}
array->SetTypedTuple(i, &tupleValues[0]);
}
array->SetName(arrayName);
grid->GetCellData()->AddArray(array);
}
//-----------------------------------------------------------------------------
void CreatePointData(vtkDataSet* grid, int numberOfComponents, int offset,
const char* arrayName)
{
vtkIdType numberOfPoints = grid->GetNumberOfPoints();
VTK_CREATE(vtkDoubleArray, array);
array->SetNumberOfComponents(numberOfComponents);
array->SetNumberOfTuples(numberOfPoints);
std::vector<double> tupleValues(numberOfComponents);
double point[3];
for(vtkIdType i=0;i<numberOfPoints;i++)
{
grid->GetPoint(i, point);
for(int j=0;j<numberOfComponents;j++)
{// +offset makes the curl/vorticity nonzero
tupleValues[j] = point[(j+offset)%3];
}
array->SetTypedTuple(i, &tupleValues[0]);
}
array->SetName(arrayName);
grid->GetPointData()->AddArray(array);
}
//-----------------------------------------------------------------------------
int IsGradientCorrect(vtkDoubleArray* gradients, int offset)
{
int numberOfComponents = gradients->GetNumberOfComponents();
for(vtkIdType i=0;i<gradients->GetNumberOfTuples();i++)
{
double* values = gradients->GetTuple(i);
for(int origComp=0;origComp<numberOfComponents/3;origComp++)
{
for(int gradDir=0;gradDir<3;gradDir++)
{
if((origComp-gradDir+offset)%3 == 0)
{
if(fabs(values[origComp*3+gradDir]-1.) > Tolerance)
{
vtkGenericWarningMacro("Gradient[ "
<< (origComp*3+gradDir)
<< " ] value should be one but is "
<< values[origComp*3+gradDir]);
if(i>10)
return 0;
}
}
else if(fabs(values[origComp*3+gradDir]) > Tolerance)
{
vtkGenericWarningMacro("Gradient[ "
<< (origComp*3+gradDir)
<< " ] value should be zero but is "
<< values[origComp*3+gradDir]);
if(i>10)
return 0;
}
}
}
}
return 1;
}
//-----------------------------------------------------------------------------
// we assume that the gradients are correct and so we can compute the "real"
// vorticity from it
int IsVorticityCorrect(vtkDoubleArray* gradients, vtkDoubleArray* vorticity)
{
if(gradients->GetNumberOfComponents() != 9 ||
vorticity->GetNumberOfComponents() != 3)
{
vtkGenericWarningMacro("Bad number of components.");
return 0;
}
for(vtkIdType i=0;i<gradients->GetNumberOfTuples();i++)
{
double* g = gradients->GetTuple(i);
double* v = vorticity->GetTuple(i);
if(!ArePointsWithinTolerance(v[0], g[7]-g[5]))
{
vtkGenericWarningMacro("Bad vorticity[0] value " << v[0] << " " <<
g[7]-g[5] << " difference is " << (v[0]-g[7]+g[5]));
return 0;
}
else if(!ArePointsWithinTolerance(v[1], g[2]-g[6]))
{
vtkGenericWarningMacro("Bad vorticity[1] value " << v[1] << " " <<
g[2]-g[6] << " difference is " << (v[1]-g[2]+g[6]));
return 0;
}
else if(!ArePointsWithinTolerance(v[2], g[3]-g[1]))
{
vtkGenericWarningMacro("Bad vorticity[2] value " << v[2] << " " <<
g[3]-g[1] << " difference is " << (v[2]-g[3]+g[1]));
return 0;
}
}
return 1;
}
//-----------------------------------------------------------------------------
// we assume that the gradients are correct and so we can compute the "real"
// Q criterion from it
int IsQCriterionCorrect(vtkDoubleArray* gradients, vtkDoubleArray* qCriterion)
{
if(gradients->GetNumberOfComponents() != 9 ||
qCriterion->GetNumberOfComponents() != 1)
{
vtkGenericWarningMacro("Bad number of components.");
return 0;
}
for(vtkIdType i=0;i<gradients->GetNumberOfTuples();i++)
{
double* g = gradients->GetTuple(i);
double qc = qCriterion->GetValue(i);
double t1 = .25*(
(g[7]-g[5])*(g[7]-g[5]) +
(g[3]-g[1])*(g[3]-g[1]) +
(g[2]-g[6])*(g[2]-g[6]) );
double t2 = .5 * ( g[0]*g[0]+g[4]*g[4]+
g[8]*g[8]+ .5 *(
(g[3]+g[1])*(g[3]+g[1]) +
(g[6]+g[2])*(g[6]+g[2]) +
(g[7]+g[5])*(g[7]+g[5]) ) );
if(!ArePointsWithinTolerance(qc, t1 - t2))
{
vtkGenericWarningMacro("Bad Q-criterion value " << qc << " " <<
t1-t2 << " difference is " << (qc-t1+t2));
return 0;
}
}
return 1;
}
//-----------------------------------------------------------------------------
// we assume that the gradients are correct and so we can compute the "real"
// divergence from it
int IsDivergenceCorrect(vtkDoubleArray* gradients, vtkDoubleArray* divergence)
{
if(gradients->GetNumberOfComponents() != 9 ||
divergence->GetNumberOfComponents() != 1)
{
vtkGenericWarningMacro("Bad number of components.");
return 0;
}
for(vtkIdType i=0;i<gradients->GetNumberOfTuples();i++)
{
double* g = gradients->GetTuple(i);
double div = divergence->GetValue(i);
double gValue = g[0]+g[4]+g[8];
if(!ArePointsWithinTolerance(div, gValue))
{
vtkGenericWarningMacro("Bad divergence value " << div << " " <<
gValue << " difference is " << (div-gValue));
return 0;
}
}
return 1;
}
//-----------------------------------------------------------------------------
int PerformTest(vtkDataSet* grid)
{
// Cleaning out the existing field data so that I can replace it with
// an analytic function that I know the gradient of
grid->GetPointData()->Initialize();
grid->GetCellData()->Initialize();
const char fieldName[] = "LinearField";
int offset = 1;
const int numberOfComponents = 3;
CreateCellData(grid, numberOfComponents, offset, fieldName);
CreatePointData(grid, numberOfComponents, offset, fieldName);
const char resultName[] = "Result";
VTK_CREATE(vtkmGradient, cellGradients);
cellGradients->SetInputData(grid);
cellGradients->SetInputScalars(
vtkDataObject::FIELD_ASSOCIATION_CELLS, fieldName);
cellGradients->SetResultArrayName(resultName);
VTK_CREATE(vtkmGradient, pointGradients);
pointGradients->SetInputData(grid);
pointGradients->SetInputScalars(
vtkDataObject::FIELD_ASSOCIATION_POINTS, fieldName);
pointGradients->SetResultArrayName(resultName);
cellGradients->Update();
pointGradients->Update();
vtkDoubleArray* gradCellArray = vtkArrayDownCast<vtkDoubleArray>(
vtkDataSet::SafeDownCast(
cellGradients->GetOutput())->GetCellData()->GetArray(resultName));
if(!grid->IsA("vtkUnstructuredGrid"))
{
// ignore cell gradients if this is an unstructured grid
// because the accuracy is so lousy
std::cout << "testing cell gradients" << std::endl;
if(!IsGradientCorrect(gradCellArray, offset))
{
return EXIT_FAILURE;
}
}
vtkDoubleArray* gradPointArray = vtkArrayDownCast<vtkDoubleArray>(
vtkDataSet::SafeDownCast(
pointGradients->GetOutput())->GetPointData()->GetArray(resultName));
std::cout << "testing point gradients" << std::endl;
if(!IsGradientCorrect(gradPointArray, offset))
{
return EXIT_FAILURE;
}
// now check on the vorticity calculations
VTK_CREATE(vtkmGradient, cellVorticity);
cellVorticity->SetInputData(grid);
cellVorticity->SetInputScalars(
vtkDataObject::FIELD_ASSOCIATION_CELLS, fieldName);
cellVorticity->SetResultArrayName(resultName);
cellVorticity->SetComputeVorticity(1);
cellVorticity->Update();
VTK_CREATE(vtkmGradient, pointVorticity);
pointVorticity->SetInputData(grid);
pointVorticity->SetInputScalars(
vtkDataObject::FIELD_ASSOCIATION_POINTS, fieldName);
pointVorticity->SetResultArrayName(resultName);
pointVorticity->SetComputeVorticity(1);
pointVorticity->SetComputeQCriterion(1);
pointVorticity->SetComputeDivergence(1);
pointVorticity->Update();
// cell stuff
vtkDoubleArray* vorticityCellArray = vtkArrayDownCast<vtkDoubleArray>(
vtkDataSet::SafeDownCast(
cellVorticity->GetOutput())->GetCellData()->GetArray("Vorticity"));
if(!IsVorticityCorrect(gradCellArray, vorticityCellArray))
{
return EXIT_FAILURE;
}
// point stuff
vtkDoubleArray* vorticityPointArray = vtkArrayDownCast<vtkDoubleArray>(
vtkDataSet::SafeDownCast(
pointVorticity->GetOutput())->GetPointData()->GetArray("Vorticity"));
if(!IsVorticityCorrect(gradPointArray, vorticityPointArray))
{
return EXIT_FAILURE;
}
vtkDoubleArray* divergencePointArray = vtkArrayDownCast<vtkDoubleArray>(
vtkDataSet::SafeDownCast(
pointVorticity->GetOutput())->GetPointData()->GetArray("Divergence"));
if(!IsDivergenceCorrect(gradPointArray, divergencePointArray))
{
return EXIT_FAILURE;
}
vtkDoubleArray* qCriterionPointArray = vtkArrayDownCast<vtkDoubleArray>(
vtkDataSet::SafeDownCast(
pointVorticity->GetOutput())->GetPointData()->GetArray("Q-criterion"));
if(!IsQCriterionCorrect(gradPointArray, qCriterionPointArray))
{
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
} // end local namespace
//-----------------------------------------------------------------------------
int TestVTKMGradientAndVorticity(int argc, char *argv[])
{
int i;
// Need to get the data root.
const char *data_root = NULL;
for (i = 0; i < argc-1; i++)