SIMSemi3D.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _SIM_SEMI3D_H
#define _SIM_SEMI3D_H
 
#include "SIMadmin.h"
#include "SIMdependency.h"
#include "SIMconfigure.h"
#include "IFEM.h"
#include "Utilities.h"
#include "MatVec.h"
#include "Vec3.h"
#include "DataExporter.h"
#include "HDF5Writer.h"
#include "tinyxml2.h"
#include <fstream>
#include <memory>
 
class TimeStep;
class VTF;
 
extern std::vector<DataExporter*> plane_exporters;
extern std::vector< std::shared_ptr<std::ostream> > plane_log_files;
 
 
template<class PlaneSolver>
class SIMSemi3D : public SIMadmin, public SIMdependency
{
public:
  typedef typename PlaneSolver::SetupProps SetupProps; 
 
  enum { dimension = 2 };
 
  explicit SIMSemi3D(const SetupProps& props_) :
    startCtx(0), planes(1), procs_per_plane(1), output_plane(-1),
    direction('Z'), props(props_)
  {
    SIMadmin::myHeading = "Plane-decoupled 3D simulation driver";
  }
 
  virtual ~SIMSemi3D()
  {
    if (!m_planes.empty())
      delete m_planes.front();
 
    for (size_t i = 1; i < m_planes.size(); i++)
    {
      // The other planes do not own the integrand, so clear the pointer to it
      // to avoid that the SIMbase destructor tries to delete it again
      m_planes[i]->clearProblem();
      delete m_planes[i];
    }
  }
 
  bool advanceStep(TimeStep& tp)
  {
    bool ok = true;
    for (size_t i = 0; i < m_planes.size() && ok; i++)
      ok = m_planes[i]->advanceStep(tp);
 
    return ok;
  }
 
  int getBDForder() const
  {
    return m_planes.empty() ? 0 : m_planes.front()->getBDForder();
  }
 
  int getDumpInterval() const
  {
    return m_planes.empty() ? 1 : m_planes.front()->getDumpInterval();
  }
 
  void printFinalNorms(const TimeStep&) {}
 
  bool preprocess(const std::vector<int>& = {}, bool = false) override
  {
    for (PlaneSolver* plane : m_planes)
      if (!plane->preprocess())
        return false;
 
    this->grabPlaneNodes();
    return true;
  }
 
  void grabPlaneNodes()
  {
    planeNodes.resize(this->getNoPlanes());
    for (size_t i = 0; i <m_planes.size(); i++)
      planeNodes[startCtx+i] = m_planes[i]->getNoNodes();
 
#ifdef HAVE_MPI
    std::vector<int> send(planeNodes);
    MPI_Allreduce(&send[0], &planeNodes[0], planeNodes.size(),
                  MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
#endif
  }
 
  std::string getName() const override { return "Semi3D"; }
 
  void registerFields(const DataExporter& exporter)
  {
    std::string name = exporter.getName();
    int plane = 1 + startCtx;
    for (size_t i = 0; i < m_planes.size(); i++, plane++)
      if (plane_exporters.size() <= i) {
        DataExporter* exp = new DataExporter(true,exporter.getStride());
        std::stringstream str;
        str << "_plane" << plane;
        HDF5Writer* hdf = new HDF5Writer(name+str.str(),
                                         m_planes[i]->getProcessAdm(),false);
        exp->registerWriter(hdf);
        plane_exporters.push_back(exp);
        m_planes[i]->registerFields(*exp);
      }
      else
        m_planes[i]->registerFields(*plane_exporters[i]);
  }
 
  bool init(const TimeStep& tp)
  {
    bool ok = true;
    for (size_t i = 0; i < m_planes.size() && ok; i++)
      ok = m_planes[i]->init(tp);
 
    return ok;
  }
 
  bool saveModel(char*,int&,int&) { return true; }
 
  bool saveStep(const TimeStep& tp, int&)
  {
    for (size_t i = 0; i < plane_exporters.size(); i++)
      plane_exporters[i]->dumpTimeLevel(&tp);
 
    // Note: VTF export is not supported for the Semi3D simulators.
    // Instead this method is used for the HDF5 export.
    return true;
  }
 
  bool serialize(std::map<std::string,std::string>&) { return false; }
  bool deSerialize(const std::map<std::string,std::string>&) { return false; }
 
  bool solveStep(TimeStep& tp)
  {
    bool ok = true;
    for (size_t i = 0; i < m_planes.size() && ok; i++) {
      m_planes[i]->getProcessAdm().cout <<"\n  Plane = "<< startCtx+i+1 <<":";
      ok = m_planes[i]->solveStep(tp);
    }
 
    return ok;
  }
 
  bool setInitialConditions()
  {
    bool ok = true;
    for (size_t i=0;i<m_planes.size();++i)
      ok &= m_planes[i]->setInitialConditions();
 
    return ok;
  }
 
  void initSystem()
  {
    for (size_t i=0;i<m_planes.size();++i)
      m_planes[i]->initSystem();
  }
 
  bool read(const char* fileName) override
  {
    if (!this->SIMadmin::read(fileName))
      return false;
 
    // Setup our communicator
#ifdef HAVE_MPI
    size_t loc_planes = planes/(nProc/procs_per_plane);
    MPI_Comm comm;
    MPI_Comm_split(PETSC_COMM_WORLD,
                   myPid/procs_per_plane,
                   myPid%procs_per_plane, &comm);
    startCtx = loc_planes*(myPid/procs_per_plane);
#else
    size_t loc_planes = planes;
#endif
 
    for (size_t i=0;i<loc_planes;++i) {
      m_planes.push_back(new PlaneSolver(props));
#ifdef HAVE_MPI
      m_planes.back()->setCommunicator(&comm);
#endif
      if (!log_files.empty()) {
        int pid = 0;
#ifdef HAVE_MPI
        MPI_Comm_rank(comm, &pid);
#endif
        if (plane_log_files.size() < loc_planes) {
          std::stringstream str;
          str << log_files <<"_plane"<< startCtx+i+1 <<".log";
          std::shared_ptr<std::ostream> file(new std::ofstream(str.str()));
          plane_log_files.push_back(file);
        }
        m_planes[i]->getProcessAdm().cout.addExtraLog(plane_log_files[i],true);
        m_planes[i]->getProcessAdm().cout.setPIDs(0, pid);
      }
      if (output_plane != -1 && output_plane != (int)(i+startCtx+1))
        m_planes[i]->getProcessAdm().cout.setNull();
      else
        m_planes[i]->getProcessAdm().cout.setStream(std::cout);
    }
 
    return true;
  }
 
  using SIMadmin::parse;
  bool parse(const tinyxml2::XMLElement* elem) override
  {
    if (!strcasecmp(elem->Value(),"postprocessing")) {
      const tinyxml2::XMLElement* child = elem->FirstChildElement();
      for (; child; child = child->NextSiblingElement())
        opt.parseOutputTag(child);
    }
    if (strcasecmp(elem->Value(),"semi3d"))
      return true;
 
    utl::getAttribute(elem, "nplanes", planes);
#ifdef HAVE_MPI
    utl::getAttribute(elem, "procs_per_plane", procs_per_plane);
    if (procs_per_plane > (size_t)nProc) procs_per_plane = nProc;
#endif
    std::string dir("Z");
    if (utl::getAttribute(elem, "direction", dir))
      direction = toupper(dir[0]);
 
    utl::getAttribute(elem,"output_prefix", log_files);
    utl::getAttribute(elem,"output_plane", output_plane);
 
    IFEM::cout <<"\tSemi3D: "<< direction
               <<" "<< planes <<" planes, "<< procs_per_plane
               <<" proces"<< (procs_per_plane > 1 ? "ses":"s") <<" per plane."
               <<"\n\tSemi3D: Printing output from ";
    if (output_plane == -1)
      IFEM::cout <<"all planes to screen."<< std::endl;
    else
      IFEM::cout <<"plane "<< output_plane <<" to screen."<< std::endl;
    if (!log_files.empty())
      IFEM::cout <<"\tSemi3D: Logging to files with prefix "
                 << log_files <<"."<< std::endl;
 
    return true;
  }
 
  using SIMdependency::registerDependency;
  template<class T>
  void registerDependency(SIMSemi3D<T>* sim, const std::string& name,
                          short int nvc, const PatchVec& patches,
                          char diffBasis = 0, int component = 0)
  {
    for (size_t i=0;i<m_planes.size(); ++i)
      if (diffBasis)
        m_planes[i]->registerDependency(sim->getPlane(i), name, nvc,
                                        sim->getPlane(i)->getFEModel(),
                                        diffBasis, component);
      else
        m_planes[i]->registerDependency(sim->getPlane(i), name, nvc);
  }
 
  size_t getNoSpaceDim() const override { return 2; }
  size_t getNoPlanes() const { return planes; }
 
  PlaneSolver* getPlane(size_t i) { return m_planes[i]; }
  const std::vector<PlaneSolver*>& getPlanes() const { return m_planes; }
 
  size_t getStartContext() const { return startCtx; }
 
  size_t getProcsPerPlane() const { return procs_per_plane; }
 
  int getPlaneNodes(int plane) const { return planeNodes[plane]; }
 
  char getDirection() const { return direction; }
 
  VTF* getVTF() { return NULL; }
  void setVTF(VTF*) {}
 
  Vec3 getForce(int) const { return Vec3(); }
  const PatchVec& getFEModel() { static PatchVec vec; return vec; }
 
  void setupDependencies()
  {
    for (size_t i=0;i<m_planes.size();++i)
      m_planes[i]->setupDependencies();
  }
 
  bool updateALE()
  {
    bool ok = true;
    for (size_t i = 0; i < m_planes.size() && ok; i++)
      ok = m_planes[i]->updateALE();
 
    return ok;
  }
 
  int getLocalNode(int) const { return -1; }
  int getGlobalNode(int) const { return -1; }
 
  int getNoBasis() const { return m_planes.front()->getNoBasis(); }
 
  size_t getNoSolutions() const { return m_planes.front()->getNoSolutions(); }
 
protected:
  std::vector<PlaneSolver*> m_planes; 
 
private:
  size_t startCtx;             
  size_t planes;               
  size_t procs_per_plane;      
  int    output_plane;         
  char   direction;            
  std::string log_files;       
  std::vector<int> planeNodes; 
  SetupProps props;            
};
 
 
template<class PlaneSolver>
struct SolverConfigurator< SIMSemi3D<PlaneSolver> > {
  int setup(SIMSemi3D<PlaneSolver>& simulator,
            const typename PlaneSolver::SetupProps& props,
            char* infile)
  {
    int retval = simulator.read(infile) ? 0 : 1;
 
    for (size_t i = 0; i < simulator.getPlanes().size() && retval == 0; i++) {
      simulator.getPlane(i)->setContext(i+simulator.getStartContext()+1);
      retval = ConfigureSIM(*simulator.getPlane(i), infile, props);
    }
    simulator.opt = simulator.getPlane(0)->opt;
 
    return retval;
  }
};
 
#endif