[cig-commits] commit 1869 by heien to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Sep 4 10:45:38 PDT 2013
Revision 1869
Simplified particle comm/IO to use only doubles, added comments for many functions/classes
U trunk/aspect/include/aspect/particle/integrator.h
U trunk/aspect/include/aspect/particle/particle.h
U trunk/aspect/include/aspect/particle/world.h
U trunk/aspect/source/particle/output.cc
U trunk/aspect/source/particle/particle.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1869&peg=1869
Diff:
Modified: trunk/aspect/include/aspect/particle/integrator.h
===================================================================
--- trunk/aspect/include/aspect/particle/integrator.h 2013-09-04 03:11:43 UTC (rev 1868)
+++ trunk/aspect/include/aspect/particle/integrator.h 2013-09-04 17:44:04 UTC (rev 1869)
@@ -29,47 +29,80 @@
{
namespace Particle
{
- // Integrator is an abstract class defining virtual methods for performing integration of
- // particle paths through the simulation velocity field
+ /**
+ * Integrator is an abstract class defining virtual methods for performing integration of
+ * particle paths through the simulation velocity field.
+ */
template <int dim, class T>
class Integrator
{
public:
virtual ~Integrator(void) {};
- // Perform an integration step of moving the particles by the specified timestep dt.
- // Implementations of this function should update the particle location.
- // If the integrator requires multiple internal steps, this should return true until
- // all internal steps are finished. Between calls to this, the particles will be moved
- // based on the velocities at their newly specified positions.
- virtual bool integrate_step(std::multimap<LevelInd, T> &particles, double dt) = 0;
+ /**
+ * Perform an integration step of moving the particles by the specified timestep dt.
+ * Implementations of this function must update the particle location.
+ * If the integrator requires multiple internal steps, this function must return
+ * true until all internal steps are finished. Between calls to this function
+ * the velocity at the updated particle positions is evaluated and passed to
+ * integrate_step during the next call.
+ *
+ * @param [in,out] particles The set of particles to integrate. The particle positions
+ * will be changed in this function based on the integration scheme.
+ * @param [in] dt The timestep length to perform the integration.
+ * @return Whether this function needs to be called again (true) for additional
+ * integration steps or if all internal steps are complete (false).
+ */
+ virtual bool integrate_step(std::multimap<LevelInd, T> &particles, const double dt) = 0;
- // Specify the MPI types and data sizes involved in transferring integration related
- // information between processes. If the integrator samples velocities at different
- // locations and the particle moves between processes during the integration step,
- // the sampled velocities must be transferred with the particle.
+ /**
+ * Secify the MPI types and data sizes involved in transferring integration related
+ * information between processes. If the integrator samples velocities at different
+ * locations and the particle moves between processes during the integration step,
+ * the sampled velocities must be transferred with the particle.
+ *
+ * @param [in,out] data_info Adds MPI data info to the specified vector indicating
+ * the quantity and type of values the integrator needs saved for this particle.
+ */
virtual void add_mpi_types(std::vector<MPIDataInfo> &data_info) = 0;
- // Must return data length in bytes of the integration related data
- // for a particle in a given format.
- virtual unsigned int data_len(ParticleDataFormat format) const = 0;
+ /**
+ * Return data length of the integration related data required for communication
+ * in terms of number of doubles.
+ *
+ * @return The number of doubles required to store the relevant integrator data.
+ */
+ virtual unsigned int data_len() const = 0;
- // Read integration related data for a particle specified by id_num
- // Returns the data pointer updated to point to the next unwritten byte
- virtual const char *read_data(ParticleDataFormat format, double id_num, const char *data) = 0;
+ /**
+ * Read integration related data for a particle specified by id_num from the data vector.
+ *
+ * @param [in] data The vector of double data to read from.
+ * @param [in] pos The position in the data vector to start reading from.
+ * @param [in] id_num The id number of the particle to read the data for.
+ * @return The position in the vector of the next unread double.
+ */
+ virtual unsigned int read_data(const std::vector<double> &data, const unsigned int &pos, const double &id_num) = 0;
- // Write integration related data for a particle specified by id_num
- // Returns the data pointer updated to point to the next unwritten byte
- virtual char *write_data(ParticleDataFormat format, double id_num, char *data) = 0;
+ /**
+ * Write integration related data to a vector for a particle
+ * specified by id_num.
+ *
+ * @param [in,out] data The vector of doubles to write integrator data into.
+ * @param [in] id_num The id number of the particle to read the data for.
+ */
+ virtual void write_data(std::vector<double> &data, const double &id_num) const = 0;
};
- // Euler scheme integrator, where y_{n+1} = y_n + dt * v(y_n)
- // This requires only one step per integration, and doesn't involve any extra data.
+ /**
+ * Euler scheme integrator, where y_{n+1} = y_n + dt * v(y_n).
+ * This requires only one step per integration, and doesn't involve any extra data.
+ */
template <int dim, class T>
class EulerIntegrator : public Integrator<dim, T>
{
public:
- virtual bool integrate_step(std::multimap<LevelInd, T> &particles, double dt)
+ virtual bool integrate_step(std::multimap<LevelInd, T> &particles, const double dt)
{
typename std::multimap<LevelInd, T>::iterator it;
Point<dim> loc, vel;
@@ -84,37 +117,38 @@
return false;
};
virtual void add_mpi_types(std::vector<MPIDataInfo> &data_info) {};
- virtual unsigned int data_len(ParticleDataFormat format) const
+ virtual unsigned int data_len() const
{
return 0;
};
- virtual const char *read_data(ParticleDataFormat format, double id_num, const char *data)
+ virtual unsigned int read_data(const std::vector<double> &data, const unsigned int &pos, const double &id_num)
{
- return data;
+ return pos;
};
- virtual char *write_data(ParticleDataFormat format, double id_num, char *data)
+ virtual void write_data(std::vector<double> &data, const double &id_num) const
{
- return data;
};
};
- // Runge Kutta second order integrator, where y_{n+1} = y_n + dt*f(k_1, k_1 = v(
- // This scheme requires storing the original location, and the read/write_data functions reflect this
+ /**
+ * Runge Kutta second order integrator, where y_{n+1} = y_n + dt*v(0.5*k_1), k_1 = dt*v(y_n).
+ * This scheme requires storing the original location, and the read/write_data functions reflect this.
+ */
template <int dim, class T>
class RK2Integrator : public Integrator<dim, T>
{
private:
- unsigned int _step;
- std::map<double, Point<dim> > _loc0;
+ unsigned int step;
+ std::map<double, Point<dim> > loc0;
public:
RK2Integrator(void)
{
- _step = 0;
- _loc0.clear();
+ step = 0;
+ loc0.clear();
};
- virtual bool integrate_step(std::multimap<LevelInd, T> &particles, double dt)
+ virtual bool integrate_step(std::multimap<LevelInd, T> &particles, const double dt)
{
typename std::multimap<LevelInd, T>::iterator it;
Point<dim> loc, vel;
@@ -125,14 +159,14 @@
id_num = it->second.get_id();
loc = it->second.get_location();
vel = it->second.get_velocity();
- if (_step == 0)
+ if (step == 0)
{
- _loc0[id_num] = loc;
+ loc0[id_num] = loc;
it->second.set_location(loc + 0.5*dt*vel);
}
- else if (_step == 1)
+ else if (step == 1)
{
- it->second.set_location(_loc0[id_num] + dt*vel);
+ it->second.set_location(loc0[id_num] + dt*vel);
}
else
{
@@ -140,93 +174,76 @@
}
}
- if (_step == 1) _loc0.clear();
- _step = (_step+1)%2;
+ if (step == 1) loc0.clear();
+ step = (step+1)%2;
// Continue until we're at the last step
- return (_step != 0);
+ return (step != 0);
};
virtual void add_mpi_types(std::vector<MPIDataInfo> &data_info)
{
- // Add the _loc0 data
- data_info.push_back(MPIDataInfo("loc0", dim, MPI_DOUBLE, sizeof(double)));
+ // Add the loc0 data
+ data_info.push_back(MPIDataInfo("loc0", dim));
};
- virtual unsigned int data_len(ParticleDataFormat format) const
+ virtual unsigned int data_len() const
{
- switch (format)
- {
- case MPI_DATA:
- case HDF5_DATA:
- return dim*sizeof(double);
- }
- return 0;
+ return dim;
};
- virtual const char *read_data(ParticleDataFormat format, double id_num, const char *data)
+ virtual unsigned int read_data(const std::vector<double> &data, const unsigned int &pos, const double &id_num)
{
- const char *p = data;
- unsigned int i;
+ unsigned int i, p = pos;
- switch (format)
+ // Read location data
+ for (i=0; i<dim; ++i)
{
- case MPI_DATA:
- case HDF5_DATA:
- // Read location data
- for (i=0; i<dim; ++i)
- {
- _loc0[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- break;
+ loc0[id_num](i) = data[p++];
}
return p;
};
- virtual char *write_data(ParticleDataFormat format, double id_num, char *data)
+ virtual void write_data(std::vector<double> &data, const double &id_num) const
{
- char *p = data;
- unsigned int i;
+ unsigned int i;
+ typename std::map<double, Point<dim> >::const_iterator it;
- // Then write our data in the appropriate format
- switch (format)
+ // Write location data
+ it = loc0.find(id_num);
+ for (i=0; i<dim; ++i)
{
- case MPI_DATA:
- case HDF5_DATA:
- // Write location data
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _loc0[id_num](i);
- p += sizeof(double);
- }
- break;
+ data.push_back(it->second(i));
}
-
- return p;
};
};
+ /**
+ * Runge Kutta fourth order integrator, where y_{n+1} = y_n + (1/6)*k1 + (1/3)*k2 + (1/3)*k3 + (1/6)*k4
+ * and k1, k2, k3, k4 are defined as usual.
+ * This scheme requires storing the original location and intermediate k1, k2, k3 values,
+ * so the read/write_data functions reflect this.
+ */
template <int dim, class T>
class RK4Integrator : public Integrator<dim, T>
{
private:
- unsigned int _step;
- std::map<double, Point<dim> > _loc0, _k1, _k2, _k3;
+ unsigned int step;
+ std::map<double, Point<dim> > loc0, k1, k2, k3;
public:
RK4Integrator(void)
{
- _step = 0;
- _loc0.clear();
- _k1.clear();
- _k2.clear();
- _k3.clear();
+ step = 0;
+ loc0.clear();
+ k1.clear();
+ k2.clear();
+ k3.clear();
};
- virtual bool integrate_step(std::multimap<LevelInd, T> &particles, double dt)
+ virtual bool integrate_step(std::multimap<LevelInd, T> &particles, const double dt)
{
typename std::multimap<LevelInd, T>::iterator it;
Point<dim> loc, vel, k4;
@@ -237,26 +254,26 @@
id_num = it->second.get_id();
loc = it->second.get_location();
vel = it->second.get_velocity();
- if (_step == 0)
+ if (step == 0)
{
- _loc0[id_num] = loc;
- _k1[id_num] = dt*vel;
- it->second.set_location(loc + 0.5*_k1[id_num]);
+ loc0[id_num] = loc;
+ k1[id_num] = dt*vel;
+ it->second.set_location(loc + 0.5*k1[id_num]);
}
- else if (_step == 1)
+ else if (step == 1)
{
- _k2[id_num] = dt*vel;
- it->second.set_location(_loc0[id_num] + 0.5*_k2[id_num]);
+ k2[id_num] = dt*vel;
+ it->second.set_location(loc0[id_num] + 0.5*k2[id_num]);
}
- else if (_step == 2)
+ else if (step == 2)
{
- _k3[id_num] = dt*vel;
- it->second.set_location(_loc0[id_num] + _k3[id_num]);
+ k3[id_num] = dt*vel;
+ it->second.set_location(loc0[id_num] + k3[id_num]);
}
- else if (_step == 3)
+ else if (step == 3)
{
k4 = dt*vel;
- it->second.set_location(_loc0[id_num] + (_k1[id_num]+2*_k2[id_num]+2*_k3[id_num]+k4)/6.0);
+ it->second.set_location(loc0[id_num] + (k1[id_num]+2*k2[id_num]+2*k3[id_num]+k4)/6.0);
}
else
{
@@ -264,117 +281,93 @@
}
}
- _step = (_step+1)%4;
- if (_step == 0)
+ step = (step+1)%4;
+ if (step == 0)
{
- _loc0.clear();
- _k1.clear();
- _k2.clear();
- _k3.clear();
+ loc0.clear();
+ k1.clear();
+ k2.clear();
+ k3.clear();
}
// Continue until we're at the last step
- return (_step != 0);
+ return (step != 0);
};
virtual void add_mpi_types(std::vector<MPIDataInfo> &data_info)
{
- // Add the _loc0, _k1, _k2, and _k3 data
- data_info.push_back(MPIDataInfo("loc0", dim, MPI_DOUBLE, sizeof(double)));
- data_info.push_back(MPIDataInfo("k1", dim, MPI_DOUBLE, sizeof(double)));
- data_info.push_back(MPIDataInfo("k2", dim, MPI_DOUBLE, sizeof(double)));
- data_info.push_back(MPIDataInfo("k3", dim, MPI_DOUBLE, sizeof(double)));
+ // Add the loc0, k1, k2, and k3 data
+ data_info.push_back(MPIDataInfo("loc0", dim));
+ data_info.push_back(MPIDataInfo("k1", dim));
+ data_info.push_back(MPIDataInfo("k2", dim));
+ data_info.push_back(MPIDataInfo("k3", dim));
};
- virtual unsigned int data_len(ParticleDataFormat format) const
+ virtual unsigned int data_len() const
{
- switch (format)
- {
- case MPI_DATA:
- case HDF5_DATA:
- return 4*dim*sizeof(double);
- }
- return 0;
+ return 4*dim;
};
- virtual const char *read_data(ParticleDataFormat format, double id_num, const char *data)
+ virtual unsigned int read_data(const std::vector<double> &data, const unsigned int &pos, const double &id_num)
{
- const char *p = data;
- unsigned int i;
+ unsigned int i, p = pos;
- switch (format)
+ // Read location data
+ for (i=0; i<dim; ++i)
{
- case MPI_DATA:
- case HDF5_DATA:
- // Read location data
- for (i=0; i<dim; ++i)
- {
- _loc0[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- // Read k1, k2 and k3
- for (i=0; i<dim; ++i)
- {
- _k1[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- _k2[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- _k3[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- break;
+ loc0[id_num](i) = data[p++];
}
+ // Read k1, k2 and k3
+ for (i=0; i<dim; ++i)
+ {
+ k1[id_num](i) = data[p++];
+ }
+ for (i=0; i<dim; ++i)
+ {
+ k2[id_num](i) = data[p++];
+ }
+ for (i=0; i<dim; ++i)
+ {
+ k3[id_num](i) = data[p++];
+ }
return p;
};
- virtual char *write_data(ParticleDataFormat format, double id_num, char *data)
+ virtual void write_data(std::vector<double> &data, const double &id_num) const
{
- char *p = data;
+ typename std::map<double, Point<dim> >::const_iterator it;
unsigned int i;
- // Then write our data in the appropriate format
- switch (format)
+ // Write location data
+ it = loc0.find(id_num);
+ for (i=0; i<dim; ++i)
{
- case MPI_DATA:
- case HDF5_DATA:
- // Write location data
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _loc0[id_num](i);
- p += sizeof(double);
- }
- // Write k1, k2 and k3
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _k1[id_num](i);
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _k2[id_num](i);
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _k3[id_num](i);
- p += sizeof(double);
- }
- break;
+ data.push_back(it->second(i));
}
-
- return p;
+ // Write k1, k2 and k3
+ it = k1.find(id_num);
+ for (i=0; i<dim; ++i)
+ {
+ data.push_back(it->second(i));
+ }
+ it = k2.find(id_num);
+ for (i=0; i<dim; ++i)
+ {
+ data.push_back(it->second(i));
+ }
+ it = k3.find(id_num);
+ for (i=0; i<dim; ++i)
+ {
+ data.push_back(it->second(i));
+ }
};
};
- // Integrator which chooses Euler, RK2 or RK4 depending on characteristics of the cell a particle is in
- // Currently used for research only
+ /**
+ * Integrator which chooses Euler, RK2 or RK4 depending on characteristics of the cell a particle is in.
+ * Currently used for research only.
+ */
template <int dim, class T>
class HybridIntegrator : public Integrator<dim, T>
{
@@ -387,13 +380,13 @@
SCHEME_RK4
};
- unsigned int _step;
- std::map<double, Point<dim> > _loc0, _k1, _k2, _k3;
- std::map<double, IntegrationScheme> _scheme;
- const parallel::distributed::Triangulation<dim> *_tria;
- const DoFHandler<dim> *_dh;
- const Mapping<dim> *_mapping;
- const TrilinosWrappers::MPI::BlockVector *_solution;
+ unsigned int step;
+ std::map<double, Point<dim> > loc0, k1, k2, k3;
+ std::map<double, IntegrationScheme> scheme;
+ const parallel::distributed::Triangulation<dim> *tria;
+ const DoFHandler<dim> *dh;
+ const Mapping<dim> *mapping;
+ const TrilinosWrappers::MPI::BlockVector *solution;
virtual IntegrationScheme select_scheme(const std::vector<Point<dim> > &cell_vertices, const std::vector<Point<dim> > &cell_velocities, const double timestep)
{
@@ -403,32 +396,31 @@
public:
HybridIntegrator(const parallel::distributed::Triangulation<dim> *new_tria, const DoFHandler<dim> *new_dh, const Mapping<dim> *new_mapping, const TrilinosWrappers::MPI::BlockVector *new_solution)
{
- _step = 0;
- _loc0.clear();
- _k1.clear();
- _k2.clear();
- _k3.clear();
- _scheme.clear();
- _tria = new_tria;
- _dh = new_dh;
- _mapping = new_mapping;
- _solution = new_solution;
+ step = 0;
+ loc0.clear();
+ k1.clear();
+ k2.clear();
+ k3.clear();
+ scheme.clear();
+ tria = new_tria;
+ dh = new_dh;
+ mapping = new_mapping;
+ solution = new_solution;
};
- virtual bool integrate_step(std::multimap<LevelInd, T> &particles, double dt)
+ virtual bool integrate_step(std::multimap<LevelInd, T> &particles, const double dt)
{
typename std::multimap<LevelInd, T>::iterator it;
Point<dim> loc, vel, k4;
double id_num;
LevelInd cur_level_ind;
IntegrationScheme cur_scheme;
- //typename parallel::distributed::Triangulation<dim>::cell_iterator found_cell;
typename DoFHandler<dim>::active_cell_iterator found_cell;
- Functions::FEFieldFunction<dim, DoFHandler<dim>, TrilinosWrappers::MPI::BlockVector> fe_value(*_dh, *_solution, *_mapping);
+ Functions::FEFieldFunction<dim, DoFHandler<dim>, TrilinosWrappers::MPI::BlockVector> fe_value(*dh, *solution, *mapping);
// If this is the first step, go through all the cells and determine
// which integration scheme the particles in each cell should use
- if (_step == 0)
+ if (step == 0)
{
Vector<double> single_res(dim+2);
std::vector<Point<dim> > cell_vertices, cell_velocities;
@@ -445,7 +437,7 @@
if (cur_level_ind != it->first)
{
cur_level_ind = it->first;
- found_cell = typename DoFHandler<dim>::active_cell_iterator(_tria, cur_level_ind.first, cur_level_ind.second, _dh);
+ found_cell = typename DoFHandler<dim>::active_cell_iterator(tria, cur_level_ind.first, cur_level_ind.second, dh);
// Ideally we should use all quadrature point velocities, but for now
// we just evaluate them at the vertices
@@ -462,7 +454,7 @@
cur_scheme = select_scheme(cell_vertices, cell_velocities, dt);
}
- _scheme[it->second.get_id()] = cur_scheme;
+ scheme[it->second.get_id()] = cur_scheme;
}
}
@@ -471,48 +463,48 @@
id_num = it->second.get_id();
loc = it->second.get_location();
vel = it->second.get_velocity();
- switch (_scheme[id_num])
+ switch (scheme[id_num])
{
case SCHEME_EULER:
- if (_step == 0)
+ if (step == 0)
{
it->second.set_location(loc + dt*vel);
}
it->second.set_vel_check(false);
break;
case SCHEME_RK2:
- if (_step == 0)
+ if (step == 0)
{
- _loc0[id_num] = loc;
+ loc0[id_num] = loc;
it->second.set_location(loc + 0.5*dt*vel);
}
- else if (_step == 1)
+ else if (step == 1)
{
- it->second.set_location(_loc0[id_num] + dt*vel);
+ it->second.set_location(loc0[id_num] + dt*vel);
}
- if (_step != 0) it->second.set_vel_check(false);
+ if (step != 0) it->second.set_vel_check(false);
break;
case SCHEME_RK4:
- if (_step == 0)
+ if (step == 0)
{
- _loc0[id_num] = loc;
- _k1[id_num] = dt*vel;
- it->second.set_location(loc + 0.5*_k1[id_num]);
+ loc0[id_num] = loc;
+ k1[id_num] = dt*vel;
+ it->second.set_location(loc + 0.5*k1[id_num]);
}
- else if (_step == 1)
+ else if (step == 1)
{
- _k2[id_num] = dt*vel;
- it->second.set_location(_loc0[id_num] + 0.5*_k2[id_num]);
+ k2[id_num] = dt*vel;
+ it->second.set_location(loc0[id_num] + 0.5*k2[id_num]);
}
- else if (_step == 2)
+ else if (step == 2)
{
- _k3[id_num] = dt*vel;
- it->second.set_location(_loc0[id_num] + _k3[id_num]);
+ k3[id_num] = dt*vel;
+ it->second.set_location(loc0[id_num] + k3[id_num]);
}
- else if (_step == 3)
+ else if (step == 3)
{
k4 = dt*vel;
- it->second.set_location(_loc0[id_num] + (_k1[id_num]+2*_k2[id_num]+2*_k3[id_num]+k4)/6.0);
+ it->second.set_location(loc0[id_num] + (k1[id_num]+2*k2[id_num]+2*k3[id_num]+k4)/6.0);
}
break;
default:
@@ -521,122 +513,97 @@
}
}
- _step = (_step+1)%4;
- if (_step == 0)
+ step = (step+1)%4;
+ if (step == 0)
{
- _loc0.clear();
- _k1.clear();
- _k2.clear();
- _k3.clear();
- _scheme.clear();
+ loc0.clear();
+ k1.clear();
+ k2.clear();
+ k3.clear();
+ scheme.clear();
}
// Continue until we're at the last step
- return (_step != 0);
+ return (step != 0);
};
virtual void add_mpi_types(std::vector<MPIDataInfo> &data_info)
{
- // Add the _loc0, _k1, _k2, and _k3 data
- data_info.push_back(MPIDataInfo("loc0", dim, MPI_DOUBLE, sizeof(double)));
- data_info.push_back(MPIDataInfo("k1", dim, MPI_DOUBLE, sizeof(double)));
- data_info.push_back(MPIDataInfo("k2", dim, MPI_DOUBLE, sizeof(double)));
- data_info.push_back(MPIDataInfo("k3", dim, MPI_DOUBLE, sizeof(double)));
- data_info.push_back(MPIDataInfo("scheme", 1, MPI_DOUBLE, sizeof(double)));
+ // Add the loc0, k1, k2, and k3 data
+ data_info.push_back(MPIDataInfo("loc0", dim));
+ data_info.push_back(MPIDataInfo("k1", dim));
+ data_info.push_back(MPIDataInfo("k2", dim));
+ data_info.push_back(MPIDataInfo("k3", dim));
+ data_info.push_back(MPIDataInfo("scheme", 1));
};
- virtual unsigned int data_len(ParticleDataFormat format) const
+ virtual unsigned int data_len() const
{
- switch (format)
- {
- case MPI_DATA:
- case HDF5_DATA:
- return (4*dim+1)*sizeof(double);
- }
- return 0;
+ return (4*dim+1);
};
- virtual const char *read_data(ParticleDataFormat format, double id_num, const char *data)
+ virtual unsigned int read_data(const std::vector<double> &data, const unsigned int &pos, const double &id_num)
{
- const char *p = data;
- unsigned int i;
+ unsigned int i, p = pos;
- switch (format)
+ // Read location data
+ for (i=0; i<dim; ++i)
{
- case MPI_DATA:
- case HDF5_DATA:
- // Read location data
- for (i=0; i<dim; ++i)
- {
- _loc0[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- // Read k1, k2 and k3
- for (i=0; i<dim; ++i)
- {
- _k1[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- _k2[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- _k3[id_num](i) = ((double *)p)[0];
- p += sizeof(double);
- }
- _scheme[id_num] = (IntegrationScheme)((double *)p)[0];
- p += sizeof(double);
- break;
+ loc0[id_num](i) = data[p++];
}
+ // Read k1, k2 and k3
+ for (i=0; i<dim; ++i)
+ {
+ k1[id_num](i) = data[p++];
+ }
+ for (i=0; i<dim; ++i)
+ {
+ k2[id_num](i) = data[p++];
+ }
+ for (i=0; i<dim; ++i)
+ {
+ k3[id_num](i) = data[p++];
+ }
+ scheme[id_num] = (IntegrationScheme)data[p++];
return p;
};
- virtual char *write_data(ParticleDataFormat format, double id_num, char *data)
+ virtual void write_data(std::vector<double> &data, const double &id_num) const
{
- char *p = data;
+ typename std::map<double, Point<dim> >::const_iterator it;
+ typename std::map<double, IntegrationScheme>::const_iterator sit;
unsigned int i;
- // Then write our data in the appropriate format
- switch (format)
+ // Write location data
+ it = loc0.find(id_num);
+ for (i=0; i<dim; ++i)
{
- case MPI_DATA:
- case HDF5_DATA:
- // Write location data
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _loc0[id_num](i);
- p += sizeof(double);
- }
- // Write k1, k2 and k3
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _k1[id_num](i);
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _k2[id_num](i);
- p += sizeof(double);
- }
- for (i=0; i<dim; ++i)
- {
- ((double *)p)[0] = _k3[id_num](i);
- p += sizeof(double);
- }
- // Write the integration scheme for this particle
- ((double *)p)[0] = _scheme[id_num];
- p += sizeof(double);
- break;
+ data.push_back(it->second(i));
}
-
- return p;
+ // Write k1, k2 and k3
+ it = k1.find(id_num);
+ for (i=0; i<dim; ++i)
+ {
+ data.push_back(it->second(i));
+ }
+ it = k2.find(id_num);
+ for (i=0; i<dim; ++i)
+ {
+ data.push_back(it->second(i));
+ }
+ it = k3.find(id_num);
+ for (i=0; i<dim; ++i)
+ {
+ data.push_back(it->second(i));
+ }
+ // Write the integration scheme for this particle
+ sit = scheme.find(id_num);
+ data.push_back((double)(sit->second));
};
};
}
}
#endif
+
Modified: trunk/aspect/include/aspect/particle/particle.h
===================================================================
--- trunk/aspect/include/aspect/particle/particle.h 2013-09-04 03:11:43 UTC (rev 1868)
+++ trunk/aspect/include/aspect/particle/particle.h 2013-09-04 17:44:04 UTC (rev 1869)
@@ -29,7 +29,9 @@
{
namespace Particle
{
- // Typedef of cell level/index pair
+ /**
+ * Typedef of cell level/index pair
+ */
typedef std::pair<int, int> LevelInd;
class MPIDataInfo
@@ -37,215 +39,296 @@
public:
std::string name;
unsigned int n_elements;
- MPI_Datatype data_type;
- unsigned int size_in_bytes;
MPIDataInfo(std::string name,
- unsigned int num_elems,
- MPI_Datatype data_type,
- unsigned int elem_size_bytes)
+ unsigned int num_elems)
:
name(name),
- n_elements(num_elems),
- data_type(data_type),
- size_in_bytes(elem_size_bytes) {};
+ n_elements(num_elems) {};
};
- enum ParticleDataFormat
- {
- MPI_DATA,
- HDF5_DATA
- };
-
- // Base class of particles - represents a particle with position, velocity, and an ID number
+ /**
+ * Base class of particles - represents a particle with position, velocity,
+ * and an ID number. This class can be extended to include data related
+ * to a particle. An example of this is shown in the DataParticle class.
+ */
template <int dim>
class BaseParticle
{
private:
- // Current particle location
+ /**
+ * Current particle location
+ */
Point<dim> location;
- // Current particle velocity
+ /**
+ * Current particle velocity
+ */
Point<dim> velocity;
- // Globally unique ID of particle
-//TODO: make this an unsigned int. but this needs adjustment in data_len, read_data, write_data and all of the places that call write_data and then parse the output...
+ /**
+ * Globally unique ID of particle
+ */
double id;
- // Whether this particle is in the local subdomain or not
+ /**
+ * Whether this particle is in the local subdomain or not
+ */
bool is_local;
- // Whether to check the velocity of this particle
- // This is used for integration schemes which require multiple
- // integration steps for some particles, but not for others
+ /**
+ * Whether to check the velocity of this particle
+ * This is used for integration schemes which require multiple
+ * integration steps for some particles, but not for others
+ */
bool check_vel;
public:
+ /**
+ * Empty constructor for BaseParticle, creates a particle at the origin with zero velocity.
+ */
BaseParticle ();
+ /**
+ * Constructor for BaseParticle, creates a particle with the specified ID at the
+ * specified location with zero velocity. Note that Aspect does not check for
+ * duplicate particle IDs so the user must be sure the IDs are unique over all processes.
+ *
+ * @param[in] new_loc Initial location of particle.
+ * @param[in] new_id Globally unique ID number of particle.
+ */
BaseParticle (const Point<dim> &new_loc,
const double &new_id);
+ /**
+ * Destructor for BaseParticle
+ */
virtual
~BaseParticle ();
+ /**
+ * Get the number of doubles required to represent this particle for communication.
+ *
+ * @return Number of doubles required to represent this particle
+ */
static unsigned int
- data_len (ParticleDataFormat format);
+ data_len ();
- virtual const char *
- read_data (ParticleDataFormat format,
- const char *data);
+ /**
+ * Read the particle data from the specified vector of doubles.
+ *
+ * @param [in] data The vector of double data to read from.
+ * @param [in] pos The position in the data vector to start reading from.
+ * @return The position in the vector of the next unread double.
+ */
+ virtual unsigned int read_data(const std::vector<double> &data, const unsigned int &pos);
+ /**
+ * Write particle data to a vector of doubles.
+ *
+ * @param [in,out] data The vector of doubles to write integrator data into.
+ */
+ virtual void write_data(std::vector<double> &data) const;
- virtual char *
- write_data (ParticleDataFormat format,
- char *data) const;
More information about the CIG-COMMITS
mailing list