Commit 4e0da8e1 authored by Burlen Loring's avatar Burlen Loring

fixes for oscillator particle miniapp

* fixes use of index space coordinates for particle in block
  tests
* support for arbitrary coordinate system (partial)
* put long functions in cpp files rather than headers
* name files by the class the define
* use diy::Grid rather than our outdated embedded copy
* remove outdated copy of diy::Grid etc
* various cleanup
parent 011a26bc
#include "Block.h"
// --------------------------------------------------------------------------
void Block::update_fields(float t)
{
// update the scalar oscilltor field
diy::for_each(grid.shape(), [&](const Vertex& v)
{
auto& gv = grid(v);
gv = 0;
auto v_global = v + Vertex(&bounds.min[0]);
for (auto& o : oscillators)
gv += o.evaluate(v_global, t);
});
// update the velocity field on the particle mesh
for (auto& particle : particles)
{
particle.velocity = { 0, 0, 0 };
for (auto& o : oscillators)
{
particle.velocity += o.evaluateGradient(particle.position, t);
}
// scale the gradient to get "units" right for velocity
particle.velocity *= velocity_scale;
}
}
// --------------------------------------------------------------------------
void Block::move_particles(float dt, const diy::Master::ProxyWithLink& cp)
{
auto link = static_cast<diy::RegularGridLink*>(cp.link());
auto particle = particles.begin();
while (particle != particles.end())
{
// update particle position
particle->position += particle->velocity * dt;
// warp position if needed
// applies periodic bci
diy::Bounds<float> wsdom = world_space_bounds(domain, origin, spacing);
for (int i = 0; i < 3; ++i)
{
if ((particle->position[i] > wsdom.max[i]) ||
(particle->position[i] < wsdom.min[i]))
{
float dm = wsdom.min[i];
float dx = wsdom.max[i] - dm;
if (fabs(dx) < 1.0e-6f)
{
particle->position[i] = dm;
}
else
{
float dp = particle->position[i] - dm;
float dpdx = dp / dx;
particle->position[i] = (dpdx - floor(dpdx))*dx + dm;
}
}
}
// check if the particle has left this block
// block bounds have ghost zones
if (!contains(domain, bounds, origin, spacing, nghost, particle->position))
{
bool enqueued = false;
// search neighbor blocks for one that now conatins this particle
for (int i = 0; i < link->size(); ++i)
{
// link bounds do not have ghost zones
if (contains(link->bounds(i), origin, spacing, particle->position))
{
/*std::cerr << "moving " << *particle << " from " << gid
<< " to " << link->target(i).gid << std::endl;*/
cp.enqueue(link->target(i), *particle);
enqueued = true;
break;
}
}
if (!enqueued)
{
std::cerr << "Error: could not find appropriate neighbor for particle: "
<< *particle << std::endl;
abort();
}
particle = particles.erase(particle);
}
else
{
++particle;
}
}
}
// --------------------------------------------------------------------------
void Block::handle_incoming_particles(const diy::Master::ProxyWithLink& cp)
{
auto link = static_cast<diy::RegularGridLink*>(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);
}
}
}
std::ostream &operator<<(std::ostream &os, const Block &b)
{
os << b.gid << ": " << b.bounds.min << " - " << b.bounds.max << std::endl;
auto it = b.particles.begin();
auto end = b.particles.end();
for (; it != end; ++it)
os << " " << *it << std::endl;
return os;
}
#ifndef BlockData_h
#define BlockData_h
#include <diy/master.hpp>
#include <diy/decomposition.hpp>
#include <diy/io/bov.hpp>
#include <diy/grid.hpp>
#include <diy/vertices.hpp>
#include "Oscillator.h"
#include "Particles.h"
#include <vector>
#include <ostream>
struct Block
{
using Vertex = diy::Grid<float,3>::Vertex;
Block(int gid_, const diy::DiscreteBounds& bounds_, const diy::DiscreteBounds& domain_,
int nghost_, const std::vector<Oscillator>& oscillators_, float velocity_scale_) :
gid(gid_),
velocity_scale(velocity_scale_),
bounds(bounds_),
domain(domain_),
nghost(nghost_),
grid(Vertex(&bounds.max[0]) - Vertex(&bounds.min[0]) + Vertex::one()),
oscillators(oscillators_)
{
// FIXME -- implement a resolution independant coordinate system. for
// now coordinate system remains tied to mesh resolution.
origin[0] = origin[1] = origin[2] = 0.0f;
spacing[0] = spacing[1] = spacing[2] = 1.0f;
}
// update scalar and vector fields
void update_fields(float t);
// update pareticle positions
void move_particles(float dt, const diy::Master::ProxyWithLink& cp);
// handle particle migration
void handle_incoming_particles(const diy::Master::ProxyWithLink& cp);
// diy memory management hooks
static void *create(){ return new Block; }
static void destroy(void* b){ delete static_cast<Block*>(b); }
int gid; // block id
float velocity_scale;
diy::DiscreteBounds bounds; // cell centered index space dimensions of this block
diy::DiscreteBounds domain; // cell centered index space dimensions of the computational domain
int nghost; // number of ghost zones
diy::Grid<float,3> grid; // container for the gridded data arrays
std::vector<Particle> particles;
std::vector<Oscillator> oscillators;
float origin[3]; // lower left most corner in the computational domain
float spacing[3]; // mesh spacing
private:
// for create; to let Master manage the blocks
Block() : gid(-1), velocity_scale(1.0f), nghost(0)
{
origin[0] = origin[1] = origin[2] = 0.0f;
spacing[0] = spacing[1] = spacing[2] = 1.0f;
}
};
// send to stream in human readbale format
std::ostream &operator<<(std::ostream &os, const Block &b);
#endif
...@@ -3,7 +3,7 @@ set(sources oscillator.cpp analysis.cpp) ...@@ -3,7 +3,7 @@ set(sources oscillator.cpp analysis.cpp)
set(libs mpi diy grid opts thread timer util) set(libs mpi diy grid opts thread timer util)
if(ENABLE_SENSEI) if(ENABLE_SENSEI)
list(APPEND sources dataadaptor.cpp bridge.cpp) list(APPEND sources bridge.cpp DataAdaptor.cpp Oscillator.cpp Particles.cpp Block.cpp)
list(APPEND libs sensei) list(APPEND libs sensei)
endif() endif()
......
#include "dataadaptor.h" #include "DataAdaptor.h"
#include "Error.h" #include "Error.h"
#include <vtkCellArray.h> #include <vtkCellArray.h>
...@@ -35,7 +35,7 @@ struct DataAdaptor::DInternals ...@@ -35,7 +35,7 @@ struct DataAdaptor::DInternals
int shape[3]; int shape[3];
int ghostLevels; int ghostLevels;
std::vector<const Particles*> ParticlesData; std::vector<const std::vector<Particle>*> ParticlesData;
vtkSmartPointer<vtkMultiBlockDataSet> ParticlesDataSet; vtkSmartPointer<vtkMultiBlockDataSet> ParticlesDataSet;
std::vector<vtkSmartPointer<vtkPolyData>> ParticleBlocks; std::vector<vtkSmartPointer<vtkPolyData>> ParticleBlocks;
}; };
...@@ -123,7 +123,7 @@ void DataAdaptor::SetBlockData(int gid, float* data) ...@@ -123,7 +123,7 @@ void DataAdaptor::SetBlockData(int gid, float* data)
} }
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
void DataAdaptor::SetParticles(int gid, const Particles& particles) void DataAdaptor::SetParticles(int gid, const std::vector<Particle> &particles)
{ {
this->Internals->ParticlesData[gid] = &particles; this->Internals->ParticlesData[gid] = &particles;
} }
...@@ -388,7 +388,7 @@ int DataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName, ...@@ -388,7 +388,7 @@ int DataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName,
continue; continue;
} }
const Particles& particles = *internals.ParticlesData[cc]; const std::vector<Particle> &particles = *internals.ParticlesData[cc];
vtkPolyData* block = internals.ParticleBlocks[cc].Get(); vtkPolyData* block = internals.ParticleBlocks[cc].Get();
auto pd = block ? block->GetPointData() : NULL; auto pd = block ? block->GetPointData() : NULL;
......
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
#include <sensei/DataAdaptor.h> #include <sensei/DataAdaptor.h>
#include "particles.h" #include "Particles.h"
class vtkDataArray; class vtkDataArray;
...@@ -33,7 +33,7 @@ public: ...@@ -33,7 +33,7 @@ public:
void SetBlockData(int gid, float* data); void SetBlockData(int gid, float* data);
/// Set particles for a specific block /// Set particles for a specific block
void SetParticles(int gid, const Particles& particles); void SetParticles(int gid, const std::vector<Particle> &particles);
// SENSEI API // SENSEI API
int GetNumberOfMeshes(unsigned int &numMeshes) override; int GetNumberOfMeshes(unsigned int &numMeshes) override;
......
#include "Oscillator.h"
#include <fstream>
#include <string>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <algorithm>
#include <diy/point.hpp>
// trim() from http://stackoverflow.com/a/217605/44738
static inline std::string &ltrim(std::string &s)
{ s.erase(s.begin(), std::find_if(s.begin(), s.end(), std::not1(std::ptr_fun<int, int>(std::isspace)))); return s; }
static inline std::string &rtrim(std::string &s)
{ s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun<int, int>(std::isspace))).base(), s.end()); return s; }
static inline std::string &trim(std::string &s) { return ltrim(rtrim(s)); }
std::vector<Oscillator> read_oscillators(std::string fn)
{
std::vector<Oscillator> res;
std::ifstream in(fn);
if (!in)
throw std::runtime_error("Unable to open " + fn);
std::string line;
while(std::getline(in, line))
{
line = trim(line);
if (line.empty() || line[0] == '#') continue;
std::istringstream iss(line);
std::string stype;
iss >> stype;
auto type = Oscillator::periodic;
if (stype == "damped")
type = Oscillator::damped;
else if (stype == "decaying")
type = Oscillator::decaying;
float x,y,z;
iss >> x >> y >> z;
float r, omega0, zeta=0.0f;
iss >> r >> omega0;
if (type == Oscillator::damped)
iss >> zeta;
res.emplace_back(Oscillator { {x,y,z}, r, omega0, zeta, type });
}
return res;
}
#ifndef oscillator_h #ifndef Oscillator_h
#define oscillator_h #define Oscillator_h
#include <fstream>
#include <string> #include <string>
#include <vector> #include <cmath>
#include <stdexcept> #include <diy/point.hpp>
#include <grid/point.h>
struct Oscillator struct Oscillator
{ {
using Vertex = grid::Point<float,3>; using Vertex = diy::Point<float,3>;
static constexpr float pi = 3.14159265358979323846; static constexpr float pi = 3.14159265358979323846;
inline float evaluate(const Vertex &v, float t) const
float evaluate(Vertex v, float t) const
{ {
t *= 2*pi; t *= 2*pi;
...@@ -27,21 +23,25 @@ struct Oscillator ...@@ -27,21 +23,25 @@ struct Oscillator
float phi = acos(zeta); float phi = acos(zeta);
float val = 1. - exp(-zeta*omega0*t) * (sin(sqrt(1-zeta*zeta)*omega0*t + phi) / sin(phi)); float val = 1. - exp(-zeta*omega0*t) * (sin(sqrt(1-zeta*zeta)*omega0*t + phi) / sin(phi));
return val * dist_damp; return val * dist_damp;
} else if (type == decaying) }
else if (type == decaying)
{ {
t += 1. / omega0; t += 1. / omega0;
float val = sin(t / omega0) / (omega0 * t); float val = sin(t / omega0) / (omega0 * t);
return val * dist_damp; return val * dist_damp;
} else if (type == periodic) }
else if (type == periodic)
{ {
t += 1. / omega0; t += 1. / omega0;
float val = sin(t / omega0); float val = sin(t / omega0);
return val * dist_damp; return val * dist_damp;
} else }
return 0; // impossible
return 0.0f; // impossible
} }
inline Vertex evaluateGradient(const Vertex& x, float t) const { Vertex evaluateGradient(const Vertex& x, float t) const
{
// let f(x, t) = this->evaluate(x,t) = o(t) * g(x) // let f(x, t) = this->evaluate(x,t) = o(t) * g(x)
// where o = oscillator, g = gaussian // where o = oscillator, g = gaussian
// therefore, df/dx = o * dg/dx // therefore, df/dx = o * dg/dx
...@@ -58,47 +58,7 @@ struct Oscillator ...@@ -58,47 +58,7 @@ struct Oscillator
enum { damped, decaying, periodic } type; enum { damped, decaying, periodic } type;
}; };
using Oscillators = std::vector<Oscillator>;
// trim() from http://stackoverflow.com/a/217605/44738
static inline std::string &ltrim(std::string &s) { s.erase(s.begin(), std::find_if(s.begin(), s.end(), std::not1(std::ptr_fun<int, int>(std::isspace)))); return s; }
static inline std::string &rtrim(std::string &s) { s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun<int, int>(std::isspace))).base(), s.end()); return s; }
static inline std::string &trim(std::string &s) { return ltrim(rtrim(s)); }
Oscillators read_oscillators(std::string fn) std::vector<Oscillator> read_oscillators(std::string fn);
{
Oscillators res;
std::ifstream in(fn);
if (!in)
throw std::runtime_error("Unable to open " + fn);
std::string line;
while(std::getline(in, line))
{
line = trim(line);
if (line.empty() || line[0] == '#') continue;
std::istringstream iss(line);
std::string stype;
iss >> stype;
auto type = Oscillator::periodic;
if (stype == "damped")
type = Oscillator::damped;
else if (stype == "decaying")
type = Oscillator::decaying;
float x,y,z;
iss >> x >> y >> z;
float r, omega0, zeta=0.0f;
iss >> r >> omega0;
if (type == Oscillator::damped)
iss >> zeta;
res.emplace_back(Oscillator { {x,y,z}, r, omega0, zeta, type });
}
return res;
}
#endif #endif
#include "Particles.h"
// --------------------------------------------------------------------------
diy::DiscreteBounds remove_ghosts(const diy::DiscreteBounds &domain,
const diy::DiscreteBounds &gbounds, int nghosts)
{
diy::DiscreteBounds bds = gbounds;
// remove ghost zones
if (nghosts > 0)
{
for (int i = 0; i < 3; ++i)
{
if (domain.min[i] != bds.min[i])
bds.min[i] += nghosts;
if (domain.max[i] != bds.max[i])
bds.max[i] -= nghosts;
}
}
return bds;
}
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
std::ostream &operator<<(std::ostream &os, const Particle &particle)
{
os << "(id " << particle.id
<< ", position: (" <<particle.position[0] << ", " << particle.position[1] << ", " << particle.position[2]
<< "), velocity: (" << particle.velocity[0] << ", " << particle.velocity[1] << ", " << particle.velocity[2]
<< "))";
return os;
}
#ifndef particles_h
#define particles_h
#include <diy/types.hpp>
#include <diy/point.hpp>
#include <random>
#include <vector>
// container for a single particle
struct Particle
{
using Vertex = diy::Point<float, 3>;
int id;
Vertex position;
Vertex velocity;
};
// put the particle in the stream in human readable format
std::ostream &operator<<(std::ostream &os, const Particle &particle);
// strips the ghost zones from the block, returns a copy
// with ghosts removed. ghost zones don't go outside of the
// computational domain.
diy::DiscreteBounds remove_ghosts(const diy::DiscreteBounds &domain,
const diy::DiscreteBounds &bounds, int nghosts);
// convert from index space into world space accounting for ghost zones
template<typename coord_type>
diy::Bounds<coord_type> world_space_bounds(const diy::DiscreteBounds &bounds,
const coord_type *origin, const coord_type *spacing)
{
// the unit of bounds member are in cell centered index space
// convert from index space bounding box into world space
diy::Bounds<coord_type> world_bounds;
world_bounds.min[0] = origin[0] + spacing[0]*bounds.min[0];
world_bounds.min[1] = origin[1] + spacing[1]*bounds.min[1];
world_bounds.min[2] = origin[2] + spacing[2]*bounds.min[2];
world_bounds.max[0] = origin[0] + spacing[0]*(bounds.max[0] + 1);
world_bounds.max[1] = origin[1] + spacing[1]*(bounds.max[1] + 1);
world_bounds.max[2] = origin[2] + spacing[2]*(bounds.max[2] + 1);
// handle planar data
for (int i = 0; i < 3; ++i)
if (bounds.min[i] == bounds.max[i])
world_bounds.max[i] = world_bounds.min[i];
return world_bounds;
}
// convert from index space into world space accounting for ghost zones
template<typename coord_type>
diy::Bounds<coord_type> world_space_bounds(const diy::DiscreteBounds &domain,
const diy::DiscreteBounds &gbounds, const coord_type *origin,
const coord_type *spacing, int nghost)
{
// remove ghost zones
diy::DiscreteBounds bounds = remove_ghosts(domain, gbounds, nghost);
// the unit of bounds member are in cell centered index space
// convert from index space bounding box into world space
return world_space_bounds(bounds, origin, spacing);
}
// gerate count particles
template<typename coord_type>
std::vector<Particle> GenerateRandomParticles(std::default_random_engine& rng,
const diy::DiscreteBounds &domain, const diy::DiscreteBounds &gbounds,
const coord_type *origin, const coord_type *spacing, int nghost,
int startId, int count)
{
diy::Bounds<coord_type> world_bounds =
world_space_bounds(domain, gbounds, origin, spacing, nghost);
std::uniform_real_distribution<coord_type> rgx(world_bounds.min[0], world_bounds.max[0]);
std::uniform_real_distribution<coord_type> rgy(world_bounds.min[1], world_bounds.max[1]);
std::uniform_real_distribution<coord_type> rgz(world_bounds.min[2], world_bounds.max[2]);
std::vector<Particle> particles;
for (int i = 0; i < count; ++i)
{
Particle p;
p.id = startId + i;
p.position = { rgx(rng), rgy(rng), rgz(rng) };
p.velocity = { 0, 0, 0 };
particles.push_back(p);
}
return particles;
}
// return true if the particle is inside this block
template<typename coord_type>
bool contains(const diy::DiscreteBounds &bounds,
const coord_type *origin, const coord_type *spacing,
const typename diy::Point<coord_type,3> &v)
{
diy::Bounds<coord_type> world_bounds =
world_space_bounds(bounds, origin, spacing);
// test if the particle is inside this block
return (v[0] >= world_bounds.min[0]) && (v[0] <= world_bounds.max[0]) &&
(v[1] >= world_bounds.min[1]) && (v[1] <= world_bounds.max[1]) &&
(v[2] >= world_bounds.min[2]) && (v[2] <= world_bounds.max[2]);
}
// return true if the particle is inside this block
template<typename coord_type>
bool contains(const diy::DiscreteBounds &domain, const diy::DiscreteBounds &gbounds,
const coord_type *origin, const coord_type *spacing, int nghost,
const typename diy::Point<coord_type,3> &v)
{
diy::Bounds<coord_type> world_bounds =
world_space_bounds(domain, gbounds, origin, spacing, nghost);
// test if the particle is inside this block
return (v[0] >= world_bounds.min[0]) && (v[0] <= world_bounds.max[0]) &&
(v[1] >= world_bounds.min[1]) && (v[1] <= world_bounds.max[1]) &&
(v[2] >= world_bounds.min[2]) && (v[2] <= world_bounds.max[2]);
}
#endif
...@@ -5,18 +5,18 @@ ...@@ -5,18 +5,18 @@
#include <diy/reduce.hpp> #include <diy/reduce.hpp>
#include <diy/partners/merge.hpp> #include <diy/partners/merge.hpp>
#include <diy/io/numpy.hpp> #include <diy/io/numpy.hpp>
#include <grid/grid.h> #include <diy/grid.hpp>
#include <grid/vertices.h> #include <diy/vertices.hpp>
#include "analysis.h" #include "analysis.h"
using GridRef = grid::GridRef<float,3>; using GridRef = diy::GridRef<float,3>;
using Vertex = GridRef::Vertex; using Vertex = GridRef::Vertex;
using Vertex4D = Vertex::UPoint; using Vertex4D = Vertex::UPoint;
struct Autocorrelation struct Autocorrelation
{ {
using Grid = grid::Grid<float,4>; using Grid = diy::Grid<float,4>;
Autocorrelation(size_t window_, int gid_, Vertex from_, Vertex to_): Autocorrelation(size_t window_, int gid_, Vertex from_, Vertex to_):
window(window_), window(window_),
...@@ -36,7 +36,7 @@ struct Autocorrelation ...@@ -36,7 +36,7 @@ struct Autocorrelation
GridRef g(data, shape); GridRef g(data, shape);
// record the values // record the values
grid::for_each(g.shape(), [&](const Vertex& v) diy::for_each(g.shape(), [&](const Vertex& v)
{ {
auto gv = g(v); auto gv = g(v);
...@@ -155,7 +155,7 @@ void analysis_final(size_t k_max, size_t nblocks) ...@@ -155,7 +155,7 @@ void analysis_final(size_t k_max, size_t nblocks)
master->foreach([](Autocorrelation* b, const diy::Master::ProxyWithLink& cp) master->foreach([](Autocorrelation* b, const diy::Master::ProxyWithLink& cp)
{ {
std::vector<float> sums(b->window, 0); std::vector<float> sums(b->window, 0);
grid::for_each(b->corr.shape(), [&](const Vertex4D& v) diy::for_each(b->corr.shape(), [&](const Vertex4D& v)
{ {
size_t w = v[3]; size_t w = v[3];
sums[w] += b->corr(v); sums[w] += b->corr(v);
...@@ -195,7 +195,7 @@ void analysis_final(size_t k_max, size_t nblocks) ...@@ -195,7 +195,7 @@ void analysis_final(size_t k_max, size_t nblocks)
MaxHeapVector maxs(b->window); MaxHeapVector maxs(b->window);
if (rp.in_link().size() == 0) if (rp.in_link().size() == 0)
{ {
grid::for_each(b->corr.shape(), [&](const Vertex4D& v) diy::for_each(b->corr.shape(), [&](const Vertex4D& v)
{