Permalink
Browse files

add mechanism to copy source data from other processors. this is one …

…of two obvious parallelization strategies, suitable for relatively small amounts of source data
  • Loading branch information...
1 parent 7ca61a1 commit 7eb6d448a90e3c9332cbb4dbb02e0aaa241ed40e @benkirk benkirk committed Dec 21, 2012
@@ -211,8 +211,8 @@ int main(int argc, char** argv)
// We now will loop over every node in the source mesh
// and add it to a source point list, along with the solution
{
- MeshBase::const_node_iterator nd = mesh_a.nodes_begin();
- MeshBase::const_node_iterator end = mesh_a.nodes_end();
+ MeshBase::const_node_iterator nd = mesh_a.local_nodes_begin();
+ MeshBase::const_node_iterator end = mesh_a.local_nodes_end();
for (; nd!=end; ++nd)
{
@@ -222,6 +222,8 @@ int main(int argc, char** argv)
}
}
+ // We have only set local values - prepare for use by gathering remote gata
+ idi.prepare_for_use();
// Create a MeshlessInterpolationFunction that uses our InverseDistanceInterpolation
// object
@@ -52,9 +52,16 @@ class MeshfreeInterpolation
public:
/**
+ * "ParallelizationStrategy"
+ */
+ enum ParallelizationStrategy { SYNC_SOURCES = 0,
+ INVALID_STRATEGY};
+ /**
* Constructor.
*/
- MeshfreeInterpolation ()
+ MeshfreeInterpolation () :
+ _is_prepared (false),
+ _parallelization_strategy (SYNC_SOURCES)
{}
/**
@@ -66,7 +73,8 @@ class MeshfreeInterpolation
/**
* Same as above, but allows you to also use stream syntax.
*/
- friend std::ostream& operator << (std::ostream& os, const MeshfreeInterpolation& mfi);
+ friend std::ostream& operator << (std::ostream& os,
+ const MeshfreeInterpolation& mfi);
/**
* Clears all internal data structures and restores to a
@@ -113,6 +121,14 @@ class MeshfreeInterpolation
const std::vector<Number> &vals);
/**
+ * Prepares data structures for use.
+ *
+ * This method is virtual so that it can be overwritten or extended as required
+ * in derived classes.
+ */
+ virtual void prepare_for_use ();
+
+ /**
* Interpolate source data at target points.
* Pure virtual, must be overriden in derived classes.
*/
@@ -121,7 +137,19 @@ class MeshfreeInterpolation
std::vector<Number> &tgt_vals) const = 0;
protected:
-
+
+ /**
+ * Gathers source points and values that have been added on other processors.
+ * Note the user is responsible for adding points only once per processor if this
+ * method is called. No attempt is made to identify duplicate points.
+ *
+ * This method is virtual so that it can be overwritten or extended as required
+ * in derived classes.
+ */
+ virtual void gather_remote_data ();
+
+ bool _is_prepared;
+ ParallelizationStrategy _parallelization_strategy;
std::vector<std::string> _names;
std::vector<Point> _src_pts;
std::vector<Number> _src_vals;
@@ -24,12 +24,15 @@
#include "libmesh/point.h"
#include "libmesh/meshfree_interpolation.h"
#include "libmesh/libmesh_logging.h"
+#include "libmesh/parallel.h"
namespace libMesh
{
+ //--------------------------------------------------------------------------------
+ // MeshfreeInterpolation methods
void MeshfreeInterpolation::print_info (std::ostream& os) const
{
os << "MeshfreeInterpolation"
@@ -50,6 +53,7 @@ namespace libMesh
void MeshfreeInterpolation::clear ()
{
+ _is_prepared = false;
_names.clear();
_src_pts.clear();
_src_vals.clear();
@@ -61,7 +65,8 @@ namespace libMesh
const std::vector<Point> &pts,
const std::vector<Number> &vals)
{
- libmesh_here();
+ _is_prepared = false;
+
libmesh_assert_equal_to (field_names.size()*pts.size(), vals.size());
// If we already have field variables, we assume we are appending.
@@ -101,7 +106,129 @@ namespace libMesh
}
+
+ void MeshfreeInterpolation::prepare_for_use ()
+ {
+ switch (_parallelization_strategy)
+ {
+ case SYNC_SOURCES:
+ this->gather_remote_data();
+ break;
+
+ case INVALID_STRATEGY:
+ libmesh_error();
+ break;
+
+ default:
+ // no preparation required for other strategies
+ break;
+ }
+ _is_prepared = true;
+ }
+
+
+
+ void MeshfreeInterpolation::gather_remote_data ()
+ {
+#ifndef LIBMESH_HAVE_MPI
+
+ // no MPI -- no-op
+ return;
+
+#else
+
+ // This function must be run on all processors at once
+ parallel_only();
+
+ START_LOG ("gather_remote_data()", "MeshfreeInterpolation");
+
+ std::vector<Real> send_buf, recv_buf;
+
+ libmesh_assert_equal_to (_src_vals.size(),
+ _src_pts.size()*this->n_field_variables());
+
+ send_buf.reserve (_src_pts.size()*(3 + this->n_field_variables()));
+
+ // Everyone packs their data at the same time
+ for (unsigned int p_idx=0, v_idx=0; p_idx<_src_pts.size(); p_idx++)
+ {
+ const Point &pt(_src_pts[p_idx]);
+
+ send_buf.push_back(pt(0));
+ send_buf.push_back(pt(1));
+ send_buf.push_back(pt(2));
+
+ for (unsigned int var=0; var<this->n_field_variables(); var++)
+ {
+ libmesh_assert_less (v_idx, _src_vals.size());
+#ifdef LIBMESH_USE_COMPLEX_NUMBERS
+ send_buf.push_back (_src_vals[v_idx].real());
+ send_buf.push_back (_src_vals[v_idx].imag());
+ v_idx++;
+
+#else
+ send_buf.push_back (_src_vals[v_idx++]);
+#endif
+ }
+ }
+
+ // Send our data to everyone else. Note that MPI-1 said you could not
+ // use the same buffer in nonblocking sends, but that restriction has
+ // recently been removed.
+ std::vector<Parallel::Request> send_request(libMesh::n_processors()-1);
+
+ // Use a tag for best practices. In debug mode parallel_only() blocks above
+ // so we can be sure there is no other shenanigarry going on, but in optimized
+ // mode there is no such guarantee - other prcoessors could be somewhere else
+ // completing some other communication, and we don't want to intercept that.
+ Parallel::MessageTag tag = Parallel::Communicator_World.get_unique_tag ( 6000 );
+
+ for (unsigned int proc=0, cnt=0; proc<libMesh::n_processors(); proc++)
+ if (proc != libMesh::processor_id())
+ CommWorld.send (proc, send_buf, send_request[cnt++], tag);
+
+ // All data has been sent. Receive remote data in any order
+ for (unsigned int comm_step=0; comm_step<(libMesh::n_processors()-1); comm_step++)
+ {
+ // blocking receive
+ CommWorld.receive (Parallel::any_source, recv_buf, tag);
+
+ // Add their data to our list
+ Point pt;
+ Number val;
+ std::vector<Real>::const_iterator it=recv_buf.begin();
+ while (it != recv_buf.end())
+ {
+ pt(0) = *it, ++it;
+ pt(1) = *it, ++it;
+ pt(2) = *it, ++it;
+
+ _src_pts.push_back(pt);
+
+ for (unsigned int var=0; var<this->n_field_variables(); var++)
+ {
+#ifdef LIBMESH_USE_COMPLEX_NUMBERS
+ val.real() = *it, ++it;
+ val.imag() = *it, ++it;
+#else
+ val = *it, ++it;
+#endif
+ _src_vals.push_back(val);
+ }
+ }
+ }
+
+ Parallel::wait (send_request);
+
+ STOP_LOG ("gather_remote_data()", "MeshfreeInterpolation");
+
+#endif // LIBMESH_HAVE_MPI
+ }
+
+
+ //--------------------------------------------------------------------------------
+ // InverseDistanceInterpolation methods
template <unsigned int KDDim>
void InverseDistanceInterpolation<KDDim>::construct_kd_tree ()
{

0 comments on commit 7eb6d44

Please sign in to comment.