LatticeMaster.hpp

Go to the documentation of this file.

//                                                         //
// Copyright (c) 2003-2013 by The University of Queensland //
// Earth Systems Science Computational Centre (ESSCC)      //
// http://www.uq.edu.au/esscc                              //
//                                                         //
// Primary Business: Brisbane, Queensland, Australia       //
// Licensed under the Open Software License version 3.0    //
// http://www.opensource.org/licenses/osl-3.0.php          //
//                                                         //

// --- Project includes ---
#include "Parallel/GetRef_cmd.h"
#include "Parallel/GeometryReader.h"

// --- STL includes ---
#include <map>
#include <set>

using namespace esys::lsm;

template <class TmplParticle>
void  CLatticeMaster::readGeometry(const std::string &fileName)
{
  console.Debug()
    << "begin CLatticeMaster::readGeometry: fileName="
    << fileName << "\n";

  // setup geometry reader and get geometry info from header file
  GeometryReader geoReader(fileName);
  GeometryInfo geoInfo = geoReader.getGeometryInfo();
        
  // check if geometry info has been set previously
  if(m_bbx_has_been_set){ 
    // check if old & new geometry are compatible
    if(!m_geo_info.isCompatible(geoInfo)){
        std::cerr << "Geometry info read from file is incompatible with previously set geometry (bounding box, circular boundaries) - Model may not run properly!" << std::endl; 
    }
    // set only the geometry version, the other geometry info remains unchanged
    m_geo_info.setLsmGeoVersion(geoInfo.getLsmGeoVersion());
  } else {
    m_geo_info=geoInfo; // set full geometry meta-info

    // set model spatial domain -> sets geometry info & initializes neighbor table in workers
    if(m_geo_info.hasAnyPeriodicDimensions()){
        setSpatialDomain(m_geo_info.getBBoxCorners()[0],m_geo_info.getBBoxCorners()[1],m_geo_info.getPeriodicDimensions());
    } else {
        setSpatialDomain(m_geo_info.getBBoxCorners()[0],m_geo_info.getBBoxCorners()[1]);
    }
    m_bbx_has_been_set=true; // flag that spatial domain has been set
  }
  
  // read particles 
  GeometryReader::ParticleIterator &particleIt = geoReader.getParticleIterator();
  console.XDebug()<< "Number of particles: " << particleIt.getNumRemaining() << "\n";
  if (geoReader.getParticleType() != "Simple") {
    throw 
      std::runtime_error(
        (std::string("Unknown particle type ") + geoReader.getParticleType()).c_str()
      );
  }

  addParticles<GeometryReader::ParticleIterator,TmplParticle>(particleIt);

  // read connections
  GeometryReader::ConnectionIterator &connectionIt =
    geoReader.getConnectionIterator();
  console.XDebug()<< "Number of connections: " << connectionIt.getNumRemaining() << "\n";

  addConnections(connectionIt);
  console.Debug()<<"end CLatticeMaster::readGeometry\n";
}

template <typename TmplVisitor>
void CLatticeMaster::visitMeshFaceReferences(const string &meshName)
{
  throw std::runtime_error("CLatticeMaster::visitMeshFaceReferences: Not implemented.");
}

template <typename TmplVisitor>
void CLatticeMaster::visitMesh2dNodeReferences(const string &meshName, TmplVisitor &visitor)
{
  console.XDebug()<<"CLatticeMaster::visitMesh2dNodeReferences( " << meshName << ")\n";
  GetNodeRefCommand cmd(getGlobalRankAndComm(),meshName);
  // broadcast command
  cmd.broadcast();
  
  // receive data (multimap)
  std::multimap<int,int> ref_mmap;
  m_tml_global_comm.gather(ref_mmap);
  // collate into set
  std::set<int> ref_set; //== this is the set of node ids == 
  for(std::multimap<int,int>::iterator iter=ref_mmap.begin();
      iter!=ref_mmap.end();
      iter++){
    ref_set.insert(iter->second);
  }
  for (std::set<int>::const_iterator it = ref_set.begin(); it != ref_set.end(); it++)
  {
    visitor.visitNodeRef(*it);
  }
  console.XDebug()<<"end CLatticeMaster::visitMesh2dNodeReferences()\n";
}

template <typename TmplVisitor>
void CLatticeMaster::visitMesh2dEdgeStress(const string &meshName, TmplVisitor &visitor)
{
  console.XDebug()<<"CLatticeMaster::visitMesh2dEdgeStress( " << meshName << ")\n";
  std::multimap<int,pair<int,Vec3> > temp_mm;
  std::map<int,Vec3> data; //=== map of id, value ===

  BroadcastCommand cmd(getGlobalRankAndComm(), CMD_GETMESH2DSTRESS);
  cmd.append(meshName.c_str());
  cmd.broadcastCommand();
  cmd.broadcastBuffer();

  // get data from slaves
  m_tml_global_comm.gather(temp_mm);
  
  // add data together
  for(std::multimap<int,pair<int,Vec3> >::iterator iter=temp_mm.begin();
      iter!=temp_mm.end();
      iter++){
    if(data.find((iter->second).first)==data.end()){ // id not in data -> insert
      data.insert(iter->second);
    } else { // id is in data -> add
      data[(iter->second).first]+=(iter->second).second;
    }
  }

  for (std::map<int,Vec3>::const_iterator it = data.begin(); it != data.end(); it++)
  {
    visitor.visitRefStressPair(it->first, it->second);
  }

  cmd.wait("visitMesh2dEdgeStress");
  console.XDebug()<<"end CLatticeMaster::visitMesh2dEdgeStress()\n";  
}

template <typename TmplVisitor>
void CLatticeMaster::visitTriMeshFaceForce(
  const string &meshName,
  TmplVisitor &visitor
)
{
  console.XDebug()<<"CLatticeMaster::visitTriMeshFaceForce( " << meshName << ")\n";
  std::multimap<int,pair<int,Vec3> > temp_mm;
  std::map<int,Vec3> data; //=== map of id, value ===

  BroadcastCommand cmd(getGlobalRankAndComm(), CMD_GETTRIMESHFORCE);
  cmd.append(meshName.c_str());
  cmd.broadcastCommand();
  cmd.broadcastBuffer();

  // get data from slaves
  m_tml_global_comm.gather(temp_mm);
  
  // add data together
  for(std::multimap<int,pair<int,Vec3> >::iterator iter=temp_mm.begin();
      iter!=temp_mm.end();
      iter++){
    if(data.find((iter->second).first)==data.end()){ // id not in data -> insert
      data.insert(iter->second);
    } else { // id is in data -> add
      data[(iter->second).first]+=(iter->second).second;
    }
  }

  for (std::map<int,Vec3>::const_iterator it = data.begin(); it != data.end(); it++)
  {
    visitor.visitRefForcePair(it->first, it->second);
  }

  cmd.wait("visitTriMeshFaceStress");
  console.XDebug()<<"end CLatticeMaster::visitTriMeshFaceForce()\n";  
}

template <typename TmplVisitor, typename TmplParticle>
void CLatticeMaster::visitParticlesOfType(
  const IdVector &particleIdVector,
  TmplVisitor &visitor
)
{
  console.Debug() << "CLatticeMaster::visitParticlesOfType: enter\n";
  typedef std::multimap<int,TmplParticle> ParticleMMap;
  ParticleMMap particleMMap;

  console.Debug()
    << "CLatticeMaster::visitParticlesOfType: broadcasting command\n";
  BroadcastCommand cmd(getGlobalRankAndComm(), CMD_GETIDPARTICLEDATA);
  cmd.broadcastCommand();

  console.Debug()
    << "CLatticeMaster::visitParticlesOfType: broadcasting particle id's\n";
  m_tml_global_comm.broadcast_cont(particleIdVector);

  console.Debug()
    << "CLatticeMaster::visitParticlesOfType:"
    << " gathering particle data from workers\n";
  m_tml_global_comm.gather_packed(particleMMap);
  console.Debug()
    << "CLatticeMaster::visitParticlesOfType:"
    << " gathered " << particleMMap.size() << " particles\n";

  console.Debug()
    << "CLatticeMaster::visitParticlesOfType:"
    << " visiting particle data\n";
  for(
    typename ParticleMMap::iterator iter=particleMMap.begin();
    iter != particleMMap.end();
    iter++
  )
  {
    iter->second.visit(visitor);
  }

  cmd.wait("visitParticles");
  console.Debug() << "CLatticeMaster::visitParticlesOfType: exit\n";
}

template <typename TmplVisitor>
void CLatticeMaster::visitParticles(
  const IdVector &particleIdVector,
  TmplVisitor &visitor
)
{
  console.Debug() << "CLatticeMaster::visitParticles: enter\n";

  if (m_particle_type == "Basic")
  {
    visitParticlesOfType<TmplVisitor,CParticle>(particleIdVector, visitor);
  }
  else if (m_particle_type == "Rot")
  {
    visitParticlesOfType<TmplVisitor,CRotParticle>(particleIdVector, visitor);
  }
  else if (m_particle_type == "RotVi")
  {
    visitParticlesOfType<TmplVisitor,CRotParticleVi>(particleIdVector, visitor);
  }
  else if (m_particle_type == "RotTherm")
  {
    visitParticlesOfType<TmplVisitor,CRotThermParticle>(particleIdVector, visitor);
  }
  else
  {
    throw
      std::runtime_error(
        std::string("Unknown particle type: ") + m_particle_type
      );
  }
  console.Debug() << "CLatticeMaster::visitParticles: exit\n";
}

template <class TmplIterator, class TmplParticle>
void CLatticeMaster::addParticles(TmplIterator &particleIt)
{
  CMPIBarrier barrier(m_global_comm);
  CMPILCmdBuffer cmd_buffer(m_global_comm,m_global_rank);

  // cast storage pointer to correct type 
  vector<TmplParticle> particleVector;
  console.XDebug()
    << "CLatticeMaster::addParticles:"
    << " Reserving vector memory for particles.\n";

  const int numBroadcastParticles = 50000;
  particleVector.reserve(numBroadcastParticles);
  console.XDebug()
    << "CLatticeMaster::addParticles:"
    << " Beginning add-particle loop..." << "\n";
  for (int i = 0; particleIt.hasNext(); i++) {
    const TmplParticle particle(particleIt.next());

    particleVector.push_back(particle);
    if (((i+1) % 5000) == 0) {
      console.XDebug()
        << "CLatticeMaster::addParticles:"
        << "Adding particle with id "
        << particleVector.rbegin()->getID() << "\n";
    }

    if (((i+1) % numBroadcastParticles) == 0)
    {
      console.XDebug() << "CLatticeMaster::addParticles:"
        << " Broadcasting receive cmd...." << "\n";
      cmd_buffer.broadcast(CMD_RECEIVEPARTICLES);
      console.XDebug()
        << "CLatticeMaster::addParticles: Broadcasting particles...." << "\n";
      m_tml_global_comm.broadcast_cont_packed(particleVector);
      barrier.wait("CLatticeMaster::addParticles: Post particle broadcast.");
      barrier.wait("CLatticeMaster::addParticles: Post Command.");
      particleVector.clear();
      particleVector.reserve(numBroadcastParticles);
    }
  }
  console.XDebug()
    << "CLatticeMaster::addParticles: Done add-particle loop..." << "\n";
  console.XDebug()
    << "CLatticeMaster::addParticles: Broadcasting final particle-receive cmd"
    << "\n";
  cmd_buffer.broadcast(CMD_RECEIVEPARTICLES);
  console.XDebug() << "CLatticeMaster::addParticles:"
    << " Broadcasting final set of particles...\n";
  m_tml_global_comm.broadcast_cont_packed(particleVector);
  barrier.wait("Post final particle broadcast.");
  barrier.wait("Post final particle-broadcast command.");
  // -- build ntable
  console.XDebug()
    << "CLatticeMaster::addParticles: "
    << "Building ntable (searchNeighbours)...\n";
  searchNeighbors(true);
  console.XDebug()
    << "CLatticeMaster::addParticles: exit\n";

}

template <class TmplIterator>
void CLatticeMaster::addConnections(TmplIterator &connectionIt)
{
  console.XDebug()
    << "CLatticeMaster::addConnections: enter\n";
  
  CMPIBarrier barrier(m_global_comm);
  CMPILCmdBuffer cmd_buffer(m_global_comm,m_global_rank);
  
  const int numBroadcastConnections = 100000;
  vector<int> connectionBuffer;
  connectionBuffer.reserve(numBroadcastConnections);
  int i = 0;
  for (; connectionIt.hasNext(); i++)
  {
    typename TmplIterator::value_type data = connectionIt.next();
    connectionBuffer.push_back(data.getTag());
    connectionBuffer.push_back(data.getP1Id());
    connectionBuffer.push_back(data.getP2Id());
    if ((i+1) % 50000 == 0)
    {
      console.XDebug() << "Adding connection number " << i << "\n";
    }
    if ((i+1) % numBroadcastConnections == 0)
    {
      console.XDebug() << "CLatticeMaster::addConnections:"
        << " Broadcasting receive cmd...." << "\n";
      cmd_buffer.broadcast(CMD_RECEIVECONNECTIONS);
      console.XDebug()
        << "CLatticeMaster::addConnections: Broadcasting connections...." << "\n";
      m_tml_global_comm.broadcast_cont_packed(connectionBuffer);
      barrier.wait("CLatticeMaster::addConnections: Post connection broadcast.");
      barrier.wait("CLatticeMaster::addConnections: Post Command.");
      connectionBuffer.clear();
      connectionBuffer.reserve(numBroadcastConnections);
    }
    //m_temp_conn[data.getTag()].push_back(data.getP1Id());
    //m_temp_conn[data.getTag()].push_back(data.getP2Id());
  }
  console.XDebug()
    << "CLatticeMaster::addConnections: Done add-connection loop..." << "\n";
  console.XDebug()
    << "CLatticeMaster::addConnections: Broadcasting final connection-receive cmd"
    << "\n";
  cmd_buffer.broadcast(CMD_RECEIVECONNECTIONS);
  console.XDebug() << "CLatticeMaster::addConnections:"
    << " Broadcasting final set of connections...\n";
  m_tml_global_comm.broadcast_cont_packed(connectionBuffer);
  barrier.wait("Post final connection broadcast.");
  barrier.wait("Post final connection-broadcast command.");

  console.XDebug()<< "Added " << i << " connections..." << "\n";
  // -- build ntable
  searchNeighbors(true);
  console.XDebug()
    << "CLatticeMaster::addConnections: exit\n";
}