Merge branch 'final-heat-diffusion-integration-squash' into 'master'

Integrate heat diffusion example

See merge request vtk/vtk-m!2195
This commit is contained in:
Valentine Peltier 2024-06-28 14:40:06 -04:00
commit 2a56651db9
8 changed files with 585 additions and 0 deletions

@ -37,6 +37,7 @@ if(VTKm_ENABLE_EXAMPLES)
add_subdirectory(cosmotools)
add_subdirectory(demo)
add_subdirectory(game_of_life)
add_subdirectory(heat_diffusion)
add_subdirectory(hello_worklet)
add_subdirectory(histogram)
add_subdirectory(ising)

@ -0,0 +1,63 @@
##============================================================================
## 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.
##============================================================================
cmake_minimum_required(VERSION 3.13)
project(Heat_Diffusion CXX)
find_package(VTKm REQUIRED)
vtkm_find_gl(OPTIONAL GL GLUT GLEW)
include(CTest)
enable_testing()
add_library(Parameters parameters.cpp)
add_library(InitalCondition InitalCondition.cpp)
target_link_libraries(InitalCondition PUBLIC vtkm_filter vtkm_io Parameters)
add_executable(HeatDiffusionExample main.cxx)
target_link_libraries(HeatDiffusionExample PUBLIC InitalCondition Parameters vtkm_filter vtkm_io)
vtkm_add_target_information(HeatDiffusionExample
DROP_UNUSED_SYMBOLS MODIFY_CUDA_FLAGS
DEVICE_SOURCES heat_diffusion.cxx Diffusion.cxx InitalCondition.cpp parameters.cpp)
if(TARGET OpenGL::GL AND
TARGET GLUT::GLUT AND
TARGET GLEW::GLEW AND
TARGET vtkm_rendering )
target_link_libraries(HeatDiffusionExample PUBLIC OpenGL::GL GLEW::GLEW GLUT::GLUT vtkm_rendering)
target_compile_definitions(HeatDiffusionExample PUBLIC -DANIMATE )
else()
MESSAGE(STATUS "Either OpenGL::GL or GLUT::GLUT or GLEW::GLEW not found: animation disabled")
endif()
##file(GLOB files "${CMAKE_CURRENT_SOURCE_DIR}/Tests/test_*.cxx")
##foreach(file ${files})
## string(REGEX REPLACE "(^.*/|\\.[^.]*$)" "" file_without_ext ${file})
## add_executable(${file_without_ext} ${file})
## target_link_libraries(${file_without_ext} PUBLIC vtkm_filter vtkm_io InitalCondition Parameters)
## vtkm_add_target_information(${file_without_ext}
## DROP_UNUSED_SYMBOLS MODIFY_CUDA_FLAGS
## DEVICE_SOURCES ${file})
## add_test(${file_without_ext}_SERIAL ${file_without_ext} -d Serial)
## add_test(${file_without_ext}_OMP ${file_without_ext} -d OpenMP)
## add_test(${file_without_ext}_TBB ${file_without_ext} -d TBB)
## add_test(${file_without_ext}_CUDA ${file_without_ext} -d Cuda)
## configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/Tests/${file_without_ext}_data.data ${CMAKE_CURRENT_BINARY_DIR}/${file_without_ext}_data.data COPYONLY)
##endforeach()
##add_test(perf_Serial HeatDiffusionExample -d Serial -p -s 20 -i 1)
##add_test(perf_OpenMP HeatDiffusionExample -d OpenMP -p -s 20 -i 1)
##add_test(perf_TBB HeatDiffusionExample -d TBB -p -s 20 -i 1)
##add_test(perf_CUDA HeatDiffusionExample -d CUDA -p -s 20 -i 1)

@ -0,0 +1,110 @@
//============================================================================
// 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.
//============================================================================
#ifndef HEAT_DIFFUSION_DIFFUSION_HPP
#define HEAT_DIFFUSION_DIFFUSION_HPP
#define NEUMMAN 0
#define DERICHLET 1
#include <vtkm/worklet/WorkletPointNeighborhood.h>
#include <vtkm/filter/FilterDataSet.h>
struct UpdateHeat : public vtkm::worklet::WorkletPointNeighborhood
{
using CountingHandle = vtkm::cont::ArrayHandleCounting<vtkm::Id>;
using ControlSignature = void(CellSetIn,
FieldInNeighborhood prevstate,
FieldIn condition,
FieldIn diffuseCoeff,
FieldOut state);
using ExecutionSignature = void(_2, _3, _4, _5);
template <typename NeighIn>
VTKM_EXEC void operator()(const NeighIn& prevstate,
const vtkm::UInt8& condition,
const vtkm::Float32& diffuseCoeff,
vtkm::Float32& state) const
{
auto current = prevstate.Get(0, 0, 0);
if (condition == NEUMMAN)
{
state = diffuseCoeff * current +
(1 - diffuseCoeff) * (0.25f * (prevstate.Get(-1, 0, 0) + prevstate.Get(0, -1, 0) +
prevstate.Get(0, 1, 0) + prevstate.Get(1, 0, 0)));
}
else if (condition == DERICHLET)
{
state = current;
}
}
};
namespace vtkm {
namespace filter {
class Diffusion : public vtkm::filter::FilterDataSet<Diffusion> {
public:
// using SupportedTypes = vtkm::List<vtkm::Float32, vtkm::UInt8, int, vtkm::Vec3f_32>;
template<typename Policy>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet &input,
vtkm::filter::PolicyBase<Policy> policy) const {
vtkm::cont::ArrayHandle<vtkm::Float32> temperature;
vtkm::cont::ArrayHandle<vtkm::Float32> prevTemperature;
vtkm::cont::ArrayHandle<vtkm::UInt8> condition;
vtkm::cont::ArrayHandle<vtkm::Float32> diffuse;
vtkm::cont::ArrayHandle<int>iteration;
const vtkm::cont::DynamicCellSet &cells = input.GetCellSet();
input.GetPointField("temperature").GetData().CopyTo(prevTemperature);
input.GetPointField("condition").GetData().CopyTo(condition);
input.GetPointField("diffuseCoeff").GetData().CopyTo(diffuse);
input.GetPointField("iteration").GetData().CopyTo(iteration);
vtkm::cont::ArrayHandle<vtkm::Float32> *ptra = &prevTemperature, *ptrb = &temperature;
//Update the temperature
for (int i = 0; i < iteration.ReadPortal().Get(0); i++)
{
this->Invoke(UpdateHeat{},
vtkm::filter::ApplyPolicyCellSet(cells, policy, *this),
*ptra,
condition,
diffuse,
*ptrb);
std::swap(ptra, ptrb);
}
vtkm::cont::DataSet output;
output.CopyStructure(input);
output.AddField(vtkm::cont::make_FieldPoint("temperature", temperature));
output.AddField(vtkm::cont::make_FieldPoint("condition", condition));
output.AddField(vtkm::cont::make_FieldPoint("diffuseCoeff", diffuse));
output.AddField(vtkm::cont::make_FieldPoint("iteration", iteration));
return output;
}
template<typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT bool DoMapField(vtkm::cont::DataSet &,
const vtkm::cont::ArrayHandle<T, StorageType> &,
const vtkm::filter::FieldMetadata &,
vtkm::filter::PolicyBase<DerivedPolicy>) { return false; }
};
}
}
#endif //HEAT_DIFFUSION_DIFFUSION_HPP

@ -0,0 +1,63 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include "InitalCondition.h"
#define NEUMMAN 0
#define DERICHLET 1
template<typename CoordType>
VTKM_EXEC void FillInitialCondition::operator()(const CoordType &coord,
vtkm::UInt8 &boundary,
vtkm::Float32 &temperature,
vtkm::Float32 &diffusion) const {
if (coord[0] == -2.0f || coord[0] == 2.0f || coord[1] == -2.0f || coord[1] == 2.0f) {
temperature = std::get<0>(parameters.temperature);
boundary = DERICHLET;
diffusion = parameters.diffuse_coeff;
} else {
float dx = (float) coord[0] - 0.0f;
float dy = (float) coord[1] - 0.0f;
float distance = dx * dx + dy * dy;
float rayon = (1.25f * 1.25f);
if (distance - rayon < 0.1f && distance - rayon > -0.1f) {
temperature = std::get<0>(parameters.temperature);
boundary = DERICHLET;
diffusion = parameters.diffuse_coeff;
} else {
temperature = std::get<1>(parameters.temperature);
boundary = NEUMMAN;
diffusion = parameters.diffuse_coeff;
}
}
}
vtkm::cont::DataSet initial_condition(const Parameters &params) {
vtkm::Id2 dimensions(params.dimension, params.dimension);
vtkm::cont::DataSet dataSet = vtkm::cont::DataSetBuilderUniform::Create(dimensions,
vtkm::Vec2f{-2.0f, -2.0f},
vtkm::Vec2f{
4.0f / (float) (dimensions[0] - 1),
4.0f /
(float) (dimensions[1] - 1)});
vtkm::cont::CoordinateSystem coords = dataSet.GetCoordinateSystem("coords");
vtkm::cont::ArrayHandle<vtkm::Float32> temperature;
vtkm::cont::ArrayHandle<vtkm::UInt8> condition;
vtkm::cont::ArrayHandle<vtkm::Float32> diffuse;
vtkm::cont::Invoker invoker;
invoker(FillInitialCondition{params}, coords, condition, temperature, diffuse);
dataSet.AddField(vtkm::cont::make_FieldPoint("temperature", temperature));
dataSet.AddField(vtkm::cont::make_FieldPoint("condition", condition));
dataSet.AddField(vtkm::cont::make_FieldPoint("diffuseCoeff", diffuse));
return dataSet;
}

@ -0,0 +1,38 @@
//============================================================================
// 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.
//============================================================================
#ifndef HEAT_DIFFUSION_INITALCONDITION_H
#define HEAT_DIFFUSION_INITALCONDITION_H
#include <utility>
#include "parameters.h"
struct FillInitialCondition : public vtkm::worklet::WorkletMapField
{
Parameters parameters;
explicit FillInitialCondition(Parameters params)
: parameters(std::move(params))
{
}
using ControlSignature = void(FieldIn, FieldOut, FieldOut, FieldOut);
using ExecutionSignature = void(_1, _2, _3, _4);
template <typename CoordType>
VTKM_EXEC void operator()(const CoordType& coord,
vtkm::UInt8& boundary,
vtkm::Float32& temperature,
vtkm::Float32& diffusion) const;
};
vtkm::cont::DataSet initial_condition(const Parameters& params);
#endif //HEAT_DIFFUSION_INITALCONDITION_H

@ -0,0 +1,198 @@
//=============================================================================
//
// 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.
//
//=============================================================================
#ifdef ANIMATE
#include <GL/glew.h>
#include <GL/glut.h>
#include <vtkm/filter/WarpScalar.h>
#include <vtkm/rendering/Actor.h>
#include <vtkm/rendering/CanvasRayTracer.h>
#include <vtkm/rendering/MapperRayTracer.h>
#include <vtkm/rendering/Scene.h>
#include <vtkm/rendering/View3D.h>
#endif
#include <vtkm/cont/DataSetFieldAdd.h>
#include <vtkm/cont/Initialize.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/cont/Timer.h>
#include "Diffusion.hpp"
#include "InitalCondition.h"
#include "parameters.h"
#include <iomanip>
#ifdef ANIMATE
void display_function(vtkm::cont::DataSet* global_ptr_data, vtkm::cont::DataSet newData, bool dim)
{
vtkm::filter::Diffusion filter;
newData = filter.Execute(*global_ptr_data);
*global_ptr_data = newData;
vtkm::rendering::CanvasRayTracer canvas(800, 800);
vtkm::rendering::MapperRayTracer mapper;
vtkm::rendering::Scene scene;
if (!dim)
{
auto normal = vtkm::cont::make_ArrayHandleConstant(vtkm::Vec3f_32{ 0.f, 0.f, 1.f },
newData.GetNumberOfPoints());
newData.AddField(vtkm::cont::make_FieldPoint("normal", normal));
vtkm::filter::WarpScalar warp(0.05f);
warp.SetUseCoordinateSystemAsField(true);
warp.SetScalarFactorField("temperature");
warp.SetFieldsToPass("temperature");
vtkm::cont::DataSet warpped = warp.Execute(newData);
vtkm::cont::ArrayHandle<vtkm::Vec3f> warpped_points;
warpped.GetPointField("warpscalar").GetData().CopyTo(warpped_points);
vtkm::rendering::Actor actor(warpped.GetCellSet(),
vtkm::cont::CoordinateSystem("warp", warpped_points),
warpped.GetPointField("temperature"),
vtkm::cont::ColorTable("Cool to warm"));
scene.AddActor(actor);
}
else
{
vtkm::rendering::Actor actor(newData.GetCellSet(),
newData.GetCoordinateSystem(),
newData.GetPointField("temperature"),
vtkm::cont::ColorTable("Cool to Warm"));
scene.AddActor(actor);
}
vtkm::rendering::View3D view(scene, mapper, canvas);
if (!dim)
{
view.GetCamera().SetLookAt(vtkm::Vec3f_32{ 0.f, 0.f, 3.0f });
view.GetCamera().Azimuth(45.0);
view.GetCamera().Elevation(40.0);
}
view.Paint();
vtkm::cont::ArrayHandle<vtkm::Vec4f_32> color_buffer = view.GetCanvas().GetColorBuffer();
const void* colorArray =
vtkm::cont::ArrayHandleBasic<vtkm::Vec4f_32>(color_buffer).GetReadPointer();
glDrawPixels(
view.GetCanvas().GetWidth(), view.GetCanvas().GetHeight(), GL_RGBA, GL_FLOAT, colorArray);
glutSwapBuffers();
}
vtkm::cont::DataSet* global_ptr_data = nullptr;
vtkm::cont::DataSet newData;
bool dim;
#endif
int main(int argc, char* argv[])
{
auto opts = vtkm::cont::InitializeOptions::DefaultAnyDevice;
vtkm::cont::InitializeResult config = vtkm::cont::Initialize(argc, argv, opts);
vtkm::cont::DataSet data;
Parameters params;
read_params(argc, argv, params);
if (params.create_matrix)
{
data = initial_condition(params);
std::cout << "Matrix size: " << params.dimension << "x" << params.dimension << std::endl;
std::cout << "Temperature outside: " << std::get<0>(params.temperature) << std::endl;
std::cout << "Temperature inside: " << std::get<1>(params.temperature) << std::endl;
std::cout << "Diffusion coefficient: " << params.diffuse_coeff << std::endl;
vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(config.Device);
}
else
{
vtkm::io::VTKDataSetReader reader(params.filename);
data = reader.ReadDataSet();
}
#ifndef ANIMATE
if (params.rendering_enable)
{
std::cout << "Animation is not available : performance is running " << std::endl;
params.rendering_enable = false;
}
#endif
if (params.rendering_enable)
{
#ifdef ANIMATE
global_ptr_data = &data;
std::vector<int> ite;
ite.push_back(10);
global_ptr_data->AddField(
vtkm::cont::make_FieldPoint("iteration", vtkm::cont::make_ArrayHandle(ite)));
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
glutInitWindowSize(800, 800);
glutCreateWindow("VTK-m Heat Diffusion");
GLenum err = glewInit();
if (GLEW_OK != err)
{
std::cout << "glewInit failed\n";
}
dim = params.twoD;
glutDisplayFunc([] { display_function(global_ptr_data, newData, dim); });
glutIdleFunc([]() { glutPostRedisplay(); });
glutMainLoop();
#endif
}
else
{
vtkm::filter::Diffusion filter;
std::vector<int> ite;
ite.push_back(params.iteration);
data.AddField(vtkm::cont::make_FieldPoint("iteration", vtkm::cont::make_ArrayHandle(ite)));
vtkm::Float32 end, start;
vtkm::cont::Timer gTimer;
gTimer.Start();
std::cout << "Number of iteration: " << params.iteration << std::endl;
start = static_cast<vtkm::Float32>(gTimer.GetElapsedTime());
data = filter.Execute(data);
end = static_cast<vtkm::Float32>(gTimer.GetElapsedTime());
std::cout << std::endl << "Execution time = " << end - start << std::endl;
double flop = ((float)params.dimension * (float)params.dimension * 8.f / (end - start)) *
(float)params.iteration;
std::cout << "MFlop = " << flop / 10e6 << std::endl;
}
return 0;
}

@ -0,0 +1,79 @@
#include "parameters.h"
#include <iostream>
#include <cstring>
void display_param()
{
std::cout << "\n\nParameters: " << std::endl;
std::cout << "\t\t-h or --help\t\t\t\t\t\tHelp\n" << std::endl;
std::cout << "\t\t-d [DEVICE]\t\t\tAny, Serial, OpenMP, TBB or Cuda\n" << std::endl;
std::cout << "\t\t-f [FILENAME]\t\t\tName of the file you want to treat\n" << std::endl;
std::cout << "\t\t-p \t\t\t\tEnable performance testing\n" << std::endl;
std::cout << "\t\t-t [TEMP_OUTSIDE] [TEMP_INSIDE]\tChange the temperature of the "
"dataset.\n\t\t\t\t\t\t\t\tDefault tempratures are 100 and 10\n"
<< std::endl;
std::cout << "\t\t-s [DIMENSION]\t\t\tChange the size of the dataset.\n\t\t\t\t\t\t\t\tDefault "
"size is 2000*2000\n"
<< std::endl;
std::cout << "\t\t-i [NB_ITERATION]\t\tChange the number of iteration for the "
"filter.\n\t\t\t\t\t\t\t\tDefault number of iteration is 100\n"
<< std::endl;
std::cout << "\t\t-c [DIFF_COEFF]\t\t\tChange the diffusion coefficient of the dataset "
".\n\t\t\t\t\t\t\t\tDefault coefficient of diffusion is 0.6\n"
<< std::endl;
std::cout << "\t\t-2d\t\t\t\t 2D rendering, default is 3D" << std::endl;
}
void read_params(int argc, char** argv, Parameters& params)
{
int i = 0;
int tmpI;
while (i < argc - 1)
{
++i;
if (!strcmp("-h", argv[i]) || !strcmp("--help", argv[i]))
{
display_param();
exit(EXIT_FAILURE);
}
if (!strcmp("-p", argv[i]))
{
params.rendering_enable = false;
continue;
}
if (!strcmp("-t", argv[i]))
{
std::get<0>(params.temperature) = strtof(argv[++i], nullptr);
std::get<1>(params.temperature) = strtof(argv[++i], nullptr);
continue;
}
if (!strcmp("-s", argv[i]))
{
tmpI = std::stoi(argv[i + 1], nullptr, 10);
if (tmpI > 9)
{
params.dimension = tmpI;
i++;
}
else
{
display_param();
exit(EXIT_FAILURE);
}
}
if (!strcmp("-c", argv[i]))
{
params.diffuse_coeff = strtof(argv[++i], nullptr);
}
if (!strcmp("-i", argv[i]))
{
params.iteration = std::stoi(argv[++i], nullptr, 10);
}
if (!strcmp("-2d", argv[i]))
{
params.twoD = true;
}
}
}

@ -0,0 +1,33 @@
//============================================================================
// 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.
//============================================================================
#ifndef PARAMETERS_H
#define PARAMETERS_H
#include <iomanip>
#include <utility>
struct Parameters
{
std::string filename;
std::tuple<float, float> temperature = std::make_tuple(100.f, 10.f);
int dimension = 2000;
int iteration = 1000;
float diffuse_coeff = 0.6f;
bool rendering_enable = true;
bool create_matrix = true;
bool twoD = false;
};
void display_param();
void read_params(int argc, char** argv, Parameters& params);
#endif