Commit 4a3897ec authored by Sujin Philip's avatar Sujin Philip Committed by Aron Helser

Add particles to oscillators

parent 6aebe1c8
......@@ -50,6 +50,12 @@ void set_data(int gid, float* data)
GlobalDataAdaptor->SetBlockData(gid, data);
}
//-----------------------------------------------------------------------------
void set_particles(int gid, const Particles& particles)
{
GlobalDataAdaptor->SetParticles(gid, particles);
}
//-----------------------------------------------------------------------------
void analyze(float time)
{
......
......@@ -4,6 +4,8 @@
#include <mpi.h>
#include <string>
#include "particles.h"
namespace bridge
{
void initialize(MPI_Comm world,
......@@ -18,6 +20,7 @@ namespace bridge
const std::string& config_file);
void set_data(int gid, float* data);
void set_particles(int gid, const Particles& particles);
void analyze(float time);
......
......@@ -10,7 +10,9 @@
#include <vtkInformation.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkObjectFactory.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkSmartPointer.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnstructuredGrid.h>
......@@ -31,6 +33,10 @@ struct DataAdaptor::DInternals
std::vector<int> DataExtent;
int shape[3];
int ghostLevels;
std::vector<const Particles*> ParticlesData;
vtkSmartPointer<vtkMultiBlockDataSet> ParticlesDataSet;
std::vector<vtkSmartPointer<vtkPolyData>> ParticleBlocks;
};
inline bool areBoundsValid(const diy::DiscreteBounds& bds)
......@@ -77,6 +83,9 @@ void DataAdaptor::Initialize(size_t nblocks, const int *shape_, int ghostLevels_
internals.shape[2] = shape_[2];
internals.ghostLevels = ghostLevels_;
internals.ParticlesData.resize(nblocks);
internals.ParticleBlocks.resize(nblocks);
this->ReleaseData();
}
......@@ -112,20 +121,26 @@ void DataAdaptor::SetBlockData(int gid, float* data)
internals.Data[gid] = data;
}
//-----------------------------------------------------------------------------
void DataAdaptor::SetParticles(int gid, const Particles& particles)
{
this->Internals->ParticlesData[gid] = &particles;
}
//-----------------------------------------------------------------------------
int DataAdaptor::GetMesh(const std::string &meshName, bool structureOnly,
vtkDataObject *&mesh)
{
if (meshName != "mesh" && meshName != "ucdmesh")
if (meshName != "mesh" && meshName != "ucdmesh" && meshName != "particles")
{
SENSEI_ERROR("the miniapp provides meshes named \"mesh\" and \"ucdmesh\""
SENSEI_ERROR("the miniapp provides meshes named \"mesh\", \"ucdmesh\", and \"particles\""
" you requested \"" << meshName << "\"")
return -1;
}
DInternals& internals = (*this->Internals);
if(meshName == "ucdmesh")
if (meshName == "ucdmesh")
{
if (!internals.uMesh)
{
......@@ -157,7 +172,7 @@ int DataAdaptor::GetMesh(const std::string &meshName, bool structureOnly,
mesh = internals.uMesh;
}
else
else if (meshName == "mesh")
{
if (!internals.Mesh)
{
......@@ -168,6 +183,20 @@ int DataAdaptor::GetMesh(const std::string &meshName, bool structureOnly,
}
mesh = internals.Mesh;
}
else if (meshName == "particles")
{
if (!internals.ParticlesDataSet)
{
internals.ParticlesDataSet = vtkSmartPointer<vtkMultiBlockDataSet>::New();
internals.ParticlesDataSet->SetNumberOfBlocks(static_cast<unsigned int>(internals.ParticlesData.size()));
for (size_t cc=0; cc < internals.ParticlesData.size(); ++cc)
{
internals.ParticlesDataSet->SetBlock(
static_cast<unsigned int>(cc), this->GetParticlesBlock(cc, structureOnly));
}
}
mesh = internals.ParticlesDataSet;
}
return 0;
}
......@@ -203,7 +232,7 @@ vtkDataObject* DataAdaptor::GetUnstructuredMesh(int gid, bool structureOnly)
}
if(structureOnly == false &&
uMesh != nullptr &&
uMesh != nullptr &&
uMesh->GetNumberOfCells() == 0)
{
// Add points.
......@@ -268,58 +297,152 @@ vtkDataObject* DataAdaptor::GetUnstructuredMesh(int gid, bool structureOnly)
return uMesh;
}
//-----------------------------------------------------------------------------
int DataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName,
int association, const std::string &arrayName)
vtkDataObject* DataAdaptor::GetParticlesBlock(int gid, bool structureOnly)
{
if ((association != vtkDataObject::FIELD_ASSOCIATION_CELLS) ||
(arrayName != "data") || (meshName != "mesh" && meshName != "ucdmesh"))
DInternals& internals = (*this->Internals);
auto particles = internals.ParticlesData[gid];
if (!particles)
{
return NULL;
}
vtkSmartPointer<vtkPolyData>& block = internals.ParticleBlocks[gid];
if (!block)
{
block = vtkSmartPointer<vtkPolyData>::New();
}
if (!structureOnly && block->GetNumberOfCells() == 0 && !particles->empty())
{
block->Initialize();
vtkNew<vtkPoints> points;
vtkNew<vtkCellArray> cells;
points->Allocate(particles->size());
cells->Allocate(particles->size());
vtkIdType pointId = 0;
for (const auto& p : *particles)
{
SENSEI_ERROR("the miniapp provides a cell centered array named \"data\" "
" on a mesh named \"mesh\"")
return -1;
points->InsertNextPoint(p.position[0], p.position[1], p.position[2]);
cells->InsertNextCell(1, &pointId);
++pointId;
}
block->SetPoints(points.Get());
block->SetVerts(cells.Get());
}
return block;
}
//-----------------------------------------------------------------------------
int DataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName,
int association, const std::string &arrayName)
{
int retVal = 1;
DInternals& internals = (*this->Internals);
vtkMultiBlockDataSet* md = vtkMultiBlockDataSet::SafeDownCast(mesh);
for (unsigned int cc=0, max=md->GetNumberOfBlocks(); cc < max; ++cc)
{
if (!internals.Data[cc]) // Exclude nullptr datasets
continue;
vtkCellData *cd = nullptr;
if(meshName == "mesh")
{
vtkSmartPointer<vtkImageData>& blockMesh = internals.BlockMesh[cc];
cd = (blockMesh? blockMesh->GetCellData() : nullptr);
}
else if(meshName == "ucdmesh")
{
vtkSmartPointer<vtkUnstructuredGrid>& uMesh = internals.UnstructuredMesh[cc];
cd = (uMesh? uMesh->GetCellData() : nullptr);
}
if (cd && !cd->GetArray(arrayName.c_str()))
if ((meshName == "mesh" || meshName == "ucdmesh") && arrayName == "data" &&
association == vtkDataObject::FIELD_ASSOCIATION_CELLS)
{
for (unsigned int cc=0, max=md->GetNumberOfBlocks(); cc < max; ++cc)
{
if (!internals.Data[cc]) // Exclude NULL datasets
{
continue;
}
vtkCellData *cd = nullptr;
if(meshName == "mesh")
{
vtkSmartPointer<vtkImageData>& blockMesh = internals.BlockMesh[cc];
cd = (blockMesh? blockMesh->GetCellData() : nullptr);
}
else if(meshName == "ucdmesh")
{
vtkSmartPointer<vtkUnstructuredGrid>& uMesh = internals.UnstructuredMesh[cc];
cd = (uMesh? uMesh->GetCellData() : nullptr);
}
const diy::DiscreteBounds &ce = internals.CellExtents[cc];
vtkIdType ncells = (ce.max[0] - ce.min[0] + 1)*
(ce.max[1] - ce.min[1] + 1)*(ce.max[2] - ce.min[2] + 1);
if (cd && !cd->GetArray(arrayName.c_str()))
{
vtkIdType ncells = (ce.max[0] - ce.min[0] + 1)*
(ce.max[1] - ce.min[1] + 1)*(ce.max[2] - ce.min[2] + 1);
vtkFloatArray *fa = vtkFloatArray::New();
fa->SetName(arrayName.c_str());
fa->SetArray(internals.Data[cc], ncells, 1);
cd->SetScalars(fa);
cd->SetActiveScalars("data");
fa->Delete();
}
}
}
else if (meshName == "particles" && association == vtkDataObject::FIELD_ASSOCIATION_POINTS &&
(arrayName == "uniqueGlobalId" || arrayName == "velocity"))
{
for (unsigned int cc=0, max=md->GetNumberOfBlocks(); cc < max; ++cc)
{
if (!internals.ParticlesData[cc])
{
continue;
}
const Particles& particles = *internals.ParticlesData[cc];
vtkFloatArray *fa = vtkFloatArray::New();
fa->SetName(arrayName.c_str());
fa->SetArray(internals.Data[cc], ncells, 1);
cd->SetScalars(fa);
cd->SetActiveScalars("data");
fa->Delete();
vtkPolyData* block = internals.ParticleBlocks[cc].Get();
auto pd = block ? block->GetPointData() : NULL;
if (pd && pd->GetArray(arrayName.c_str()) == NULL)
{
if (arrayName == "uniqueGlobalId")
{
vtkNew<vtkIdTypeArray> array;
array->SetName(arrayName.c_str());
array->SetNumberOfComponents(1);
array->SetNumberOfTuples(block->GetNumberOfPoints());
for (vtkIdType i = 0; i < array->GetNumberOfTuples(); ++i)
{
vtkIdType ugid = particles[i].id;
array->SetTypedTuple(i, &ugid);
}
pd->AddArray(array.Get());
}
else if (arrayName == "velocity")
{
vtkNew<vtkFloatArray> array;
array->SetName(arrayName.c_str());
array->SetNumberOfComponents(3);
array->SetNumberOfTuples(block->GetNumberOfPoints());
float vel[3] = { 0, 0, 0 };
for (vtkIdType i = 0; i < array->GetNumberOfTuples(); ++i)
{
vel[0] = particles[i].velocity[0];
vel[1] = particles[i].velocity[1];
vel[2] = particles[i].velocity[2];
array->SetTypedTuple(i, vel);
}
pd->SetScalars(array.Get());
pd->SetActiveScalars(array->GetName());
}
}
retVal = pd ? 0 : 1;
}
}
#ifndef NDEBUG
else
{
SENSEI_ERROR("the miniapp provides a cell centered array named \"data\" "
"on meshes named \"mesh\" or \"ucdmesh\"; or point centered arrays named "
"\"uniqueGlobalId\" and \"velocity\" on a mesh named \"particles\"")
return 1;
}
#endif
return 0;
return retVal;
}
//----------------------------------------------------------------------------
int DataAdaptor::GetMeshHasGhostCells(const std::string &/*meshName*/,
int DataAdaptor::GetMeshHasGhostCells(const std::string &/*meshName*/,
int &nLayers)
{
DInternals& internals = (*this->Internals);
......@@ -446,7 +569,7 @@ int DataAdaptor::AddGhostCellsArray(vtkDataObject *mesh, const std::string &mesh
//-----------------------------------------------------------------------------
int DataAdaptor::GetNumberOfMeshes(unsigned int &numMeshes)
{
numMeshes = 2;
numMeshes = 3;
return 0;
}
......@@ -463,6 +586,11 @@ int DataAdaptor::GetMeshName(unsigned int id, std::string &meshName)
meshName = "ucdmesh";
return 0;
}
else if (id == 2)
{
meshName = "particles";
return 0;
}
SENSEI_ERROR("Failed to get mesh name")
return -1;
......@@ -478,6 +606,11 @@ int DataAdaptor::GetNumberOfArrays(const std::string &meshName, int association,
numberOfArrays = 1;
return 0;
}
else if ((meshName == "particles") && (association == vtkDataObject::POINT))
{
numberOfArrays = 2;
return 0;
}
return 0;
}
......@@ -491,6 +624,12 @@ int DataAdaptor::GetArrayName(const std::string &meshName, int association,
arrayName = "data";
return 0;
}
else if ((meshName == "particles") && (association == vtkDataObject::POINT))
{
static const char *names[] = { "uniqueGlobalId", "velocity" };
arrayName = names[index];
return 0;
}
SENSEI_ERROR("Failed to get array name")
return -1;
......@@ -500,8 +639,14 @@ int DataAdaptor::GetArrayName(const std::string &meshName, int association,
int DataAdaptor::ReleaseData()
{
DInternals& internals = (*this->Internals);
<<<<<<< HEAD
internals.Mesh = nullptr;
internals.uMesh = nullptr;
=======
internals.Mesh = NULL;
internals.uMesh = NULL;
internals.ParticlesDataSet = NULL;
>>>>>>> Add particles to oscillators
for (auto i : internals.CellExtents)
{
i.min[0] = i.min[1] = i.min[2] = 0;
......@@ -509,9 +654,17 @@ int DataAdaptor::ReleaseData()
}
for (size_t cc=0, max = internals.Data.size(); cc < max; ++cc)
{
<<<<<<< HEAD
internals.Data[cc] = nullptr;
internals.BlockMesh[cc] = nullptr;
internals.UnstructuredMesh[cc] = nullptr;
=======
internals.Data[cc] = NULL;
internals.BlockMesh[cc] = NULL;
internals.UnstructuredMesh[cc] = NULL;
internals.ParticlesData[cc] = NULL;
internals.ParticleBlocks[cc] = NULL;
>>>>>>> Add particles to oscillators
}
return 0;
}
......
......@@ -3,6 +3,8 @@
#include <sensei/DataAdaptor.h>
#include "particles.h"
class vtkDataArray;
namespace oscillators
......@@ -30,6 +32,9 @@ public:
/// Set data for a specific block.
void SetBlockData(int gid, float* data);
/// Set particles for a specific block
void SetParticles(int gid, const Particles& particles);
// SENSEI API
int GetNumberOfMeshes(unsigned int &numMeshes) override;
......@@ -61,6 +66,8 @@ protected:
vtkDataObject* GetUnstructuredMesh(int gid, bool structureOnly);
vtkDataArray* CreateGhostCellsArray(int cc) const;
vtkDataObject* GetParticlesBlock(int gid, bool structureOnly);
private:
DataAdaptor(const DataAdaptor&); // not implemented.
void operator=(const DataAdaptor&); // not implemented.
......
#include <vector>
#include <chrono>
#include <ctime>
#include <format.h>
#include <opts/opts.h>
......@@ -12,6 +13,7 @@
#include <grid/vertices.h>
#include "oscillator.h"
#include "particles.h"
#include "senseiConfig.h"
#ifdef ENABLE_SENSEI
......@@ -31,15 +33,22 @@ using Link = diy::RegularGridLink;
using Master = diy::Master;
using Proxy = Master::ProxyWithLink;
using RandomSeedType = std::default_random_engine::result_type;
inline bool IsVertexInsideBounds(const Vertex& v, const Bounds& b)
{
return (v[0] >= static_cast<float>(b.min[0])) && (v[0] <= static_cast<float>(b.max[0])) &&
(v[1] >= static_cast<float>(b.min[1])) && (v[1] <= static_cast<float>(b.max[1])) &&
(v[2] >= static_cast<float>(b.min[2])) && (v[2] <= static_cast<float>(b.max[2]));
}
struct Block
{
Block(int gid_, const Bounds& bounds_, const Vertex& domain_shape_, const Oscillators& oscillators_):
Block(int gid_, const Bounds& bounds_, const Bounds& domain_, const Oscillators* oscillators_):
gid(gid_),
bounds(bounds_),
domain_shape(domain_shape_),
from(bounds.min),
domain(domain_),
grid(Vertex(bounds.max) - Vertex(bounds.min) + Vertex::one()),
oscillators(oscillators_)
{
......@@ -51,17 +60,98 @@ struct Block
{
auto& gv = grid(v);
gv = 0;
auto v_global = v + from;
auto v_global = v + bounds.min;
for (auto& o : oscillators)
gv += o.evaluate(v_global, t);
});
for (auto& particle : particles)
{
particle.velocity = { 0, 0, 0 };
if (oscillators)
{
for (auto& o : *oscillators)
{
particle.velocity += o.evaluateGradient(particle.position, t);
}
}
// normalize
particle.velocity /= std::sqrt(particle.velocity.norm());
}
}
void move_particles(float dt, const Proxy& cp)
{
auto link = static_cast<Link*>(cp.link());
auto particle = particles.begin();
while (particle != particles.end())
{
particle->position += particle->velocity * dt;
// warp position if needed
for (int i = 0; i < 3; ++i)
{
if (particle->position[i] > static_cast<float>(domain.max[i]))
{
particle->position[i] += static_cast<float>(domain.min[i] - domain.max[i]);
}
else if (particle->position[i] < static_cast<float>(domain.min[i]))
{
particle->position[i] += static_cast<float>(domain.max[i] - domain.min[i]);
}
}
if (!IsVertexInsideBounds(particle->position, bounds))
{
bool enqueued = false;
for (int i = 0; i < link->size(); ++i)
{
if (IsVertexInsideBounds(particle->position, link->bounds(i)))
{
cp.enqueue(link->target(i), *particle);
enqueued = true;
break;
}
}
if (!enqueued)
{
fmt::print(
"Error: could not find appropriate neighbor for particle: "
"id: {}, position: ({}, {}, {}), velocity: ({}, {}, {})\n",
particle->id,
particle->position[0], particle->position[1], particle->position[2],
particle->velocity[0], particle->velocity[1], particle->velocity[2]);
}
particle = particles.erase(particle);
}
else
{
++particle;
}
}
}
void handle_incoming_particles(const Proxy& cp)
{
auto link = static_cast<Link*>(cp.link());
for (int i = 0; i < link->size(); ++i)
{
auto nbr = link->target(i).gid;
while(cp.incoming(nbr))
{
Particle particle;
cp.dequeue(nbr, particle);
particles.push_back(particle);
}
}
}
void analyze_block()
{
#ifdef ENABLE_SENSEI
bridge::set_data(gid, grid.data());
bridge::set_particles(gid, particles);
#else
analyze(gid, grid.data());
#endif
......@@ -72,10 +162,10 @@ struct Block
int gid;
Bounds bounds;
Vertex domain_shape;
Vertex from;
Bounds domain;
Grid grid;
Oscillators oscillators;
Particles particles;
const Oscillators* oscillators;
private:
Block() {} // for create; to let Master manage the blocks
......@@ -97,6 +187,8 @@ int main(int argc, char** argv)
size_t k_max = 3;
int threads = 1;
int ghostLevels = 0;
int numberOfParticles = 1024;
int seed = -1;
std::string config_file;
std::string out_prefix = "";
Options ops(argc, argv);
......@@ -113,7 +205,9 @@ int main(int argc, char** argv)
>> Option( "t-end", t_end, "end time")
>> Option('j', "jobs", threads, "number of threads to use")
>> Option('o', "output", out_prefix, "prefix to save output")
>> Option('g', "ghost levels", ghostLevels, "Number of ghost levels")
>> Option('g', "ghost levels", ghostLevels, "number of ghost levels")
>> Option('p', "particles", numberOfParticles, "number of random particles to generate")
>> Option( "seed", seed, "specify a random seed")
;
bool sync = ops >> Present("sync", "synchronize after each time step");
bool log = ops >> Present("log", "generate full time and memory usage log");
......@@ -142,6 +236,21 @@ int main(int argc, char** argv)
return 1;
}
int particlesPerBlock = (numberOfParticles + nblocks - 1) / nblocks;
if (seed == -1)
{
if (world.rank() == 0)
{
seed = static_cast<int>(std::time(nullptr));
fmt::print("using seed {}\n", seed);
}
diy::mpi::broadcast(world, seed, 0);
}