From 81a90a0b41a1316f2411f297587b215d4a8ec4a4 Mon Sep 17 00:00:00 2001 From: Dave Pugmire Date: Fri, 21 Sep 2018 15:02:23 -0400 Subject: [PATCH 1/9] Add shell code for particle advection filter. --- src/vtkh/filters/CMakeLists.txt | 8 +- src/vtkh/filters/ParticleAdvection.cpp | 110 +++++++++++++++++++++++++ src/vtkh/filters/ParticleAdvection.hpp | 39 +++++++++ 3 files changed, 154 insertions(+), 3 deletions(-) create mode 100644 src/vtkh/filters/ParticleAdvection.cpp create mode 100644 src/vtkh/filters/ParticleAdvection.hpp diff --git a/src/vtkh/filters/CMakeLists.txt b/src/vtkh/filters/CMakeLists.txt index 4a42d615..9dc6b3db 100644 --- a/src/vtkh/filters/CMakeLists.txt +++ b/src/vtkh/filters/CMakeLists.txt @@ -12,7 +12,8 @@ set(vtkh_filters_headers NoOp.hpp Lagrangian.hpp MarchingCubes.hpp - PointAverage.hpp + ParticleAdvection.hpp + PointAverage.hpp Recenter.hpp Threshold.hpp Slice.hpp @@ -29,6 +30,7 @@ set(vtkh_filters_sources NoOp.cpp Lagrangian.cpp MarchingCubes.cpp + ParticleAdvection.cpp PointAverage.cpp Recenter.cpp Threshold.cpp @@ -80,14 +82,14 @@ if (MPI_FOUND) # triggers cuda compile list(APPEND vtkh_filters_mpi_deps cuda) endif() - + blt_add_library( NAME vtkh_filters_mpi SOURCES ${vtkh_filters_sources} HEADERS ${vtkh_filters_headers} DEPENDS_ON ${vtkh_filters_mpi_deps} ) - + if(ENABLE_OPENMP) if(CUDA_FOUND) blt_add_target_compile_flags(TO vtkh_filters_mpi FLAGS "-Xcompiler ${OpenMP_CXX_FLAGS}") diff --git a/src/vtkh/filters/ParticleAdvection.cpp b/src/vtkh/filters/ParticleAdvection.cpp new file mode 100644 index 00000000..94f4988a --- /dev/null +++ b/src/vtkh/filters/ParticleAdvection.cpp @@ -0,0 +1,110 @@ +#include +#include +#include +#include +#include + +namespace vtkh +{ + +ParticleAdvection::ParticleAdvection() +{ + +} + +ParticleAdvection::~ParticleAdvection() +{ + +} + +void +ParticleAdvection::SetField(const std::string &field_name) +{ + m_field_name = field_name; +} + +void +ParticleAdvection::SetStepSize(const double &step_size) +{ + m_step_size = step_size; +} + +void +ParticleAdvection::SetWriteFrequency(const int &write_frequency) +{ + m_write_frequency = write_frequency; +} + +void +ParticleAdvection::SetCustomSeedResolution(const int &cust_res) +{ + m_cust_res = cust_res; +} + +void +ParticleAdvection::SetSeedResolutionInX(const int &x_res) +{ + m_x_res = x_res; +} + +void +ParticleAdvection::SetSeedResolutionInY(const int &y_res) +{ + m_y_res = y_res; +} +void +ParticleAdvection::SetSeedResolutionInZ(const int &z_res) +{ + m_z_res = z_res; +} + + +void ParticleAdvection::PreExecute() +{ + Filter::PreExecute(); +} + +void ParticleAdvection::PostExecute() +{ + Filter::PostExecute(); +} + +void ParticleAdvection::DoExecute() +{ + vtkm::worklet::ParticleAdvection particleAdvectionWorklet; + + this->m_output = new DataSet(); + const int num_domains = this->m_input->GetNumberOfDomains(); + + for(int i = 0; i < num_domains; ++i) + { + vtkm::Id domain_id; + vtkm::cont::DataSet dom; + this->m_input->GetDomain(i, dom, domain_id); + if(dom.HasField(m_field_name)) + { + using vectorField_d = vtkm::cont::ArrayHandle>; + using vectorField_f = vtkm::cont::ArrayHandle>; + auto field = dom.GetField(m_field_name).GetData(); + if(!field.IsSameType(vectorField_d()) && !field.IsSameType(vectorField_f())) + { + throw Error("Vector field type does not match > or >"); + } + } + else + { + throw Error("Domain does not contain specified vector field for ParticleAdvection analysis."); + } + + vtkm::cont::DataSet res; + m_output->AddDomain(res, domain_id); + } +} + +std::string +ParticleAdvection::GetName() const +{ + return "vtkh::ParticleAdvection"; +} + +} // namespace vtkh diff --git a/src/vtkh/filters/ParticleAdvection.hpp b/src/vtkh/filters/ParticleAdvection.hpp new file mode 100644 index 00000000..fa8e248a --- /dev/null +++ b/src/vtkh/filters/ParticleAdvection.hpp @@ -0,0 +1,39 @@ +#ifndef VTK_H_PARTICLE_ADVECTION_HPP +#define VTK_H_PARTICLE_ADVECTION_HPP + +#include +#include +#include + +namespace vtkh +{ + +class ParticleAdvection : public Filter +{ +public: + ParticleAdvection(); + virtual ~ParticleAdvection(); + std::string GetName() const override; + void SetField(const std::string &field_name); + void SetStepSize(const double &step_size); + void SetWriteFrequency(const int &write_frequency); + void SetCustomSeedResolution(const int &cust_res); + void SetSeedResolutionInX(const int &x_res); + void SetSeedResolutionInY(const int &y_res); + void SetSeedResolutionInZ(const int &z_res); + + +protected: + void PreExecute() override; + void PostExecute() override; + void DoExecute() override; + + std::string m_field_name; + double m_step_size; + int m_write_frequency; + int m_cust_res; + int m_x_res, m_y_res, m_z_res; +}; + +} //namespace vtkh +#endif From 4132adf37c1cd48b47596ef198e2b4a600e6bdfc Mon Sep 17 00:00:00 2001 From: Dave Pugmire Date: Wed, 21 Nov 2018 13:03:18 -0500 Subject: [PATCH 2/9] Initial code for doing distributed memory streamlines. --- src/vtkh/DataSet.cpp | 120 ++--- src/vtkh/filters/BoundsMap.hpp | 70 +++ src/vtkh/filters/CMakeLists.txt | 14 +- src/vtkh/filters/CommData.hpp | 64 +++ src/vtkh/filters/Communicator.hpp | 253 ++++++++++ src/vtkh/filters/DebugMeowMeow.hpp | 13 + src/vtkh/filters/Integrator.hpp | 241 ++++++++++ src/vtkh/filters/Lagrangian.cpp | 43 +- src/vtkh/filters/MemStream.h | 239 ++++++++++ src/vtkh/filters/MemStream.hxx | 518 ++++++++++++++++++++ src/vtkh/filters/Particle.hpp | 36 ++ src/vtkh/filters/ParticleAdvection.cpp | 360 ++++++++++++-- src/vtkh/filters/ParticleAdvection.hpp | 99 +++- src/vtkh/filters/adapter.h | 23 + src/vtkh/filters/avtParICAlgorithm.hpp | 114 +++++ src/vtkh/filters/avtParICAlgorithm.hxx | 636 +++++++++++++++++++++++++ src/vtkh/filters/util.hpp | 129 +++++ 17 files changed, 2833 insertions(+), 139 deletions(-) create mode 100644 src/vtkh/filters/BoundsMap.hpp create mode 100644 src/vtkh/filters/CommData.hpp create mode 100644 src/vtkh/filters/Communicator.hpp create mode 100644 src/vtkh/filters/DebugMeowMeow.hpp create mode 100644 src/vtkh/filters/Integrator.hpp create mode 100644 src/vtkh/filters/MemStream.h create mode 100644 src/vtkh/filters/MemStream.hxx create mode 100644 src/vtkh/filters/Particle.hpp create mode 100644 src/vtkh/filters/adapter.h create mode 100644 src/vtkh/filters/avtParICAlgorithm.hpp create mode 100644 src/vtkh/filters/avtParICAlgorithm.hxx create mode 100644 src/vtkh/filters/util.hpp diff --git a/src/vtkh/DataSet.cpp b/src/vtkh/DataSet.cpp index 7aa80150..69736eb9 100644 --- a/src/vtkh/DataSet.cpp +++ b/src/vtkh/DataSet.cpp @@ -23,7 +23,7 @@ bool GlobalAgreement(bool local) { bool agreement = local; #ifdef VTKH_PARALLEL - int local_boolean = local ? 1 : 0; + int local_boolean = local ? 1 : 0; int global_boolean; MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle()); MPI_Allreduce((void *)(&local_boolean), @@ -55,13 +55,13 @@ class MemSetWorklet : public vtkm::worklet::WorkletMapField typedef void ControlSignature(FieldOut<>); typedef void ExecutionSignature(_1); - + VTKM_EXEC void operator()(T &value) const { value = Value; } -}; //class MemSetWorklet +}; //class MemSetWorklet template void MemSet(vtkm::cont::ArrayHandle &array, const T value, const vtkm::Id num_values) @@ -73,8 +73,8 @@ void MemSet(vtkm::cont::ArrayHandle &array, const T value, const vtkm::Id num } // namespace detail -void -DataSet::AddDomain(vtkm::cont::DataSet data_set, vtkm::Id domain_id) +void +DataSet::AddDomain(vtkm::cont::DataSet data_set, vtkm::Id domain_id) { if(m_domains.size() != 0) { @@ -87,7 +87,7 @@ DataSet::AddDomain(vtkm::cont::DataSet data_set, vtkm::Id domain_id) m_domain_ids.push_back(domain_id); } -vtkm::cont::Field +vtkm::cont::Field DataSet::GetField(const std::string &field_name, const vtkm::Id domain_index) { assert(domain_index >= 0); @@ -97,7 +97,7 @@ DataSet::GetField(const std::string &field_name, const vtkm::Id domain_index) } vtkm::cont::DataSet& -DataSet::GetDomain(const vtkm::Id index) +DataSet::GetDomain(const vtkm::Id index) { const size_t num_domains = m_domains.size(); @@ -108,7 +108,7 @@ DataSet::GetDomain(const vtkm::Id index) <<" in "<(m_domains.size()); } -vtkm::Id +vtkm::Id DataSet::GetGlobalNumberOfDomains() const { - vtkm::Id domains = this->GetNumberOfDomains(); -#ifdef VKTH_PARALLEL + vtkm::Id domains = this->GetNumberOfDomains(); +#ifdef VTKH_PARALLEL MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle()); - int local_doms = static_cast(domains); + int local_doms = static_cast(domains); int global_doms = 0; - MPI_Allreduce(&local_doms, - &global_doms, - 1, - MPI_INT, + MPI_Allreduce(&local_doms, + &global_doms, + 1, + MPI_INT, MPI_SUM, mpi_comm); domains = global_doms; @@ -164,7 +164,7 @@ DataSet::GetGlobalNumberOfDomains() const return domains; } -vtkm::Bounds +vtkm::Bounds DataSet::GetDomainBounds(const int &domain_index, vtkm::Id coordinate_system_index) const { @@ -172,8 +172,8 @@ DataSet::GetDomainBounds(const int &domain_index, vtkm::cont::CoordinateSystem coords; try { - coords = m_domains[domain_index].GetCoordinateSystem(index); - } + coords = m_domains[domain_index].GetCoordinateSystem(index); + } catch (const vtkm::cont::Error &error) { std::stringstream msg; @@ -187,7 +187,7 @@ DataSet::GetDomainBounds(const int &domain_index, } -vtkm::Bounds +vtkm::Bounds DataSet::GetBounds(vtkm::Id coordinate_system_index) const { const vtkm::Id index = coordinate_system_index; @@ -204,7 +204,7 @@ DataSet::GetBounds(vtkm::Id coordinate_system_index) const return bounds; } -vtkm::Bounds +vtkm::Bounds DataSet::GetGlobalBounds(vtkm::Id coordinate_system_index) const { vtkm::Bounds bounds; @@ -227,7 +227,7 @@ DataSet::GetGlobalBounds(vtkm::Id coordinate_system_index) const vtkm::Float64 global_z_max = 0; MPI_Allreduce((void *)(&x_min), - (void *)(&global_x_min), + (void *)(&global_x_min), 1, MPI_DOUBLE, MPI_MIN, @@ -241,7 +241,7 @@ DataSet::GetGlobalBounds(vtkm::Id coordinate_system_index) const mpi_comm); MPI_Allreduce((void *)(&y_min), - (void *)(&global_y_min), + (void *)(&global_y_min), 1, MPI_DOUBLE, MPI_MIN, @@ -255,7 +255,7 @@ DataSet::GetGlobalBounds(vtkm::Id coordinate_system_index) const mpi_comm); MPI_Allreduce((void *)(&z_min), - (void *)(&global_z_min), + (void *)(&global_z_min), 1, MPI_DOUBLE, MPI_MIN, @@ -278,7 +278,7 @@ DataSet::GetGlobalBounds(vtkm::Id coordinate_system_index) const return bounds; } -vtkm::cont::ArrayHandle +vtkm::cont::ArrayHandle DataSet::GetRange(const std::string &field_name) const { const size_t num_domains = m_domains.size(); @@ -296,12 +296,12 @@ DataSet::GetRange(const std::string &field_name) const const vtkm::cont::Field &field = m_domains[i].GetField(field_name); vtkm::cont::ArrayHandle sub_range; sub_range = field.GetRange(); - - num_components = sub_range.GetPortalConstControl().GetNumberOfValues(); + + num_components = sub_range.GetPortalConstControl().GetNumberOfValues(); range = sub_range; - vtkm::Id components = sub_range.GetPortalConstControl().GetNumberOfValues(); - + vtkm::Id components = sub_range.GetPortalConstControl().GetNumberOfValues(); + // first range with data. Set range and keep looking if(num_components == 0) { @@ -309,7 +309,7 @@ DataSet::GetRange(const std::string &field_name) const range = sub_range; continue; } - + // This is not the first valid range encountered. // Validate and expand the current range if(components != num_components) @@ -329,11 +329,11 @@ DataSet::GetRange(const std::string &field_name) const range.GetPortalControl().Set(c, c_range); } } - + return range; } -vtkm::cont::ArrayHandle +vtkm::cont::ArrayHandle DataSet::GetGlobalRange(const std::string &field_name) const { vtkm::cont::ArrayHandle range; @@ -346,7 +346,7 @@ DataSet::GetGlobalRange(const std::string &field_name) const // it is possible to have an empty dataset at one of the ranks // so we must check for this so MPI comm does not hang. // We also want to check for num components mis-match - // + // int *global_components = new int[vtkh::GetMPISize()]; int comps = static_cast(num_components); @@ -369,7 +369,7 @@ DataSet::GetGlobalRange(const std::string &field_name) const components = global_components[i]; continue; } - + // verify that this matches are current components if(global_components[i] != 0 && components != global_components[i]) { @@ -404,12 +404,12 @@ DataSet::GetGlobalRange(const std::string &field_name) const local_min = std::numeric_limits::max(); local_max = std::numeric_limits::min(); } - + vtkm::Float64 global_min = 0; vtkm::Float64 global_max = 0; MPI_Allreduce((void *)(&local_min), - (void *)(&global_min), + (void *)(&global_min), 1, MPI_DOUBLE, MPI_MIN, @@ -432,7 +432,7 @@ DataSet::GetGlobalRange(const std::string &field_name) const return range; } -void +void DataSet::PrintSummary(std::ostream &stream) const { for(size_t dom = 0; dom < m_domains.size(); ++dom) @@ -442,7 +442,7 @@ DataSet::PrintSummary(std::ostream &stream) const } } -bool +bool DataSet::IsPointMesh(const vtkm::Id cell_set_index) const { bool is_points = true; @@ -459,7 +459,7 @@ DataSet::IsPointMesh(const vtkm::Id cell_set_index) const return is_points; } -bool +bool DataSet::IsStructured(int &topological_dims, const vtkm::Id cell_set_index) const { topological_dims = -1; @@ -468,14 +468,14 @@ DataSet::IsStructured(int &topological_dims, const vtkm::Id cell_set_index) cons for(size_t i = 0; i < num_domains; ++i) { const vtkm::cont::DataSet &dom = m_domains[i]; - int dims; + int dims; is_structured = VTKMDataSetInfo::IsStructured(dom, dims, cell_set_index) && is_structured; if(i == 0) { - topological_dims = dims; + topological_dims = dims; } - + if(!is_structured || dims != topological_dims) { topological_dims = -1; @@ -501,7 +501,7 @@ DataSet::SetCycle(const vtkm::UInt64 cycle) vtkm::UInt64 DataSet::GetCycle() const { - return m_cycle; + return m_cycle; } DataSet::DataSet() @@ -513,8 +513,8 @@ DataSet::~DataSet() { } -vtkm::cont::DataSet& -DataSet::GetDomainById(const vtkm::Id domain_id) +vtkm::cont::DataSet& +DataSet::GetDomainById(const vtkm::Id domain_id) { const size_t size = m_domain_ids.size(); @@ -536,11 +536,11 @@ bool DataSet::HasDomainId(const vtkm::Id &domain_id) const { if(m_domain_ids[i] == domain_id) return true; } - + return false; } -void +void DataSet::AddConstantPointField(const vtkm::Float32 value, const std::string fieldname) { const size_t size = m_domain_ids.size(); @@ -555,7 +555,7 @@ DataSet::AddConstantPointField(const vtkm::Float32 value, const std::string fiel } } -bool +bool DataSet::FieldExists(const std::string &field_name) const { bool exists = false; @@ -572,12 +572,12 @@ DataSet::FieldExists(const std::string &field_name) const return exists; } -bool +bool DataSet::GlobalFieldExists(const std::string &field_name) const { bool exists = FieldExists(field_name); #ifdef VTKH_VTKH_PARALLEL - int local_boolean = exists ? 1 : 0; + int local_boolean = exists ? 1 : 0; int global_boolean; MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle()); @@ -608,11 +608,11 @@ DataSet::GetFieldAssociation(const std::string field_name, bool &valid_field) co valid_field = true; if(!this->GlobalFieldExists(field_name)) { - valid_field = false; + valid_field = false; return vtkm::cont::Field::Association::ANY; } - - int assoc_id = -1; + + int assoc_id = -1; if(this->FieldExists(field_name)) { const size_t num_domains = m_domains.size(); @@ -649,7 +649,7 @@ DataSet::GetFieldAssociation(const std::string field_name, bool &valid_field) co } } } - + #ifdef VTKH_PARALLEL MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle()); diff --git a/src/vtkh/filters/BoundsMap.hpp b/src/vtkh/filters/BoundsMap.hpp new file mode 100644 index 00000000..baf80f45 --- /dev/null +++ b/src/vtkh/filters/BoundsMap.hpp @@ -0,0 +1,70 @@ +#ifndef VTK_H_BOUNDS_MAP_HPP +#define VTK_H_BOUNDS_MAP_HPP + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +class BoundsMap +{ +public: + BoundsMap() {} + BoundsMap(const BoundsMap &_bm) : bm(_bm.bm) {} + + void AddBlock(int id, const vtkm::Bounds &bounds) + { + if (bm.find(id) == bm.end()) + bm[id] = bounds; + else + throw "Duplicate block"; + } + template