mainappmpi.cpp

Go to the documentation of this file.

#include "mainappmpi.h"
#include "unittests.h"
#include "hdf5orthowriter.h"
#include <stdarg.h>
#include <string.h>
#include <libgen.h>
#include "netmpi.h"
#include "fcompoundvector.h"
#include "timer.h"
#include "dummyflash.h"
/*
#include "laxfriedrichssystemmpi.h"
#include "hybridpressionmoduleagm_mpi.h"
#include "hdf5orthowriter.h"
#include "russanovsystemmpi.h"
#include "hybridpressionmodulempi.h"
#include "laxfriedrichsmethodmpi.h"
#include "russanovmethodmpi.h"
#include "hybridpressionmoduleumfpack_mpi.h"
#include "dummyflash.h"
*/


MainAppMPI::MainAppMPI() 
   :MainApp()
{
  if (NetMPI::nProcess() == 1)
    setOutput(std::cout);
  else
  {
    char strOut[200];
    if (getRank() == 0)
      system("rm out-*");
    NetMPI::Barrier();
    sprintf(strOut,"out-%d",getRank());
    out.open(strOut);
    setOutput(out);
    NetMPI::setLog(out);
  }
}





 MainAppMPI::~MainAppMPI() 
{

}



void MainAppMPI::setLocalTriangulation() 
{

  if (this->runningInParallel())
  {
    if (getLocalNElems()[0] < 2*configFile.getUnsigned("MESH_OVERLAP_SIZE"))
    {
      throw new Exception("Mesh Overlap gets the entire domain");
    }
  }
  Point<3> P0 = getLocalP0();
  Point<3> P1 = getLocalP1();

  std::vector<unsigned> v = getLocalNElems();
  VecWellInfo wells = getWells();

  
  //OK generate the grid
  m_pMesh = new OrthoMesh(P0,P1,v[0],v[1],v[2]);
  m_pMesh->putWells(wells);

  
}


int MainAppMPI::getNProcess() 
{
  int nP;
  MPI_Comm_size(MPI_COMM_WORLD,&nP);
  return nP;
}

int MainAppMPI::getRank() 
{
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
  return rank;
}


Point3D MainAppMPI::getLocalP0() 
{
  std::vector<unsigned> nElems = getNElements();
  Point3D DX = getDX();

  int rank = getRank();
  int nProcess = getNProcess();

  //Get the first startElement of the non overlapping region
  unsigned startElement = nElems[0]/nProcess*rank;
  int rest = nElems[0]%nProcess;
  startElement +=( (rank <  rest) ?  rank : rest);

  //Now add the overlap region
  if (rank != 0)
  {
    startElement-=configFile.getInt("MESH_OVERLAP_SIZE");
  }
  return Point3D(startElement*DX(0),0,0);
}


Point3D MainAppMPI::getLocalP1() 
{
  Point3D LP1  = MainApp::getP1();
  Point3D DX = getDX();
  LP1(0) = getLocalP0()[0]+getLocalNElems()[0]*DX(0);
  return LP1;
}



std::vector<unsigned> MainAppMPI::getLocalNElems() 
{
  std::vector<unsigned> v(3);

  int rank     = getRank();
  int nProcess = getNProcess();
  std::vector<unsigned> nElems = MainApp::getNElements();


  int rest = nElems[0]%nProcess;
  int nElements =  nElems[0]/nProcess + (rank < rest ? 1 : 0);
  
  if (nProcess == 1)
    nElements+=0;
  else if (rank == 0 || rank == (nProcess -1))
    nElements+=configFile.getInt("MESH_OVERLAP_SIZE");
  else  
    nElements+=2*configFile.getInt("MESH_OVERLAP_SIZE");


    
  v[0]=nElements;
  v[1]=nElems[1];
  v[2]=nElems[2];
  return v;
}


void MainAppMPI::execute(std::string fileName) 
{
  try {
    //Read the config file.
    configFile.readFile(fileName);
        
    if (configFile.getInt("DEBUG_PRINT_TRIANGULATION",0) == 1)
    {
      getMesh().print();
    }

    if (configFile.getInt("UNIT_TESTS",0))
    {
      UnitTests tests(*this);
      tests.execute();
      printBeginSection("END UNITS TESTS");
      return;
    }
    
    //Initialize HDF5 writer with the triangulation.
    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
    hdf5.setOutputFile(getHDF5FileName());
    print("HDF5","HDF5 Ouput: %s",getHDF5FileName().c_str());
    hdf5.writeMesh("tria0",getMesh());
    hdf5.setVariable("Time",0);
    
    m_pSequencer = getSequencer();
      
    switch (configFile.getInt("SEQUENCE_ALGORITHM"))
    {
    case 0:
      print("Sequencer","Dynamic Module Test\n");
      printLog(std::cout);
      m_pTransportModule = getTransportModule();
      m_pDynamicModule = getDynamicModule();
      m_pSequencer->testDynamicModule(*m_pDynamicModule,*m_pTransportModule);
      break;
    case 1:
      {
        m_pTransportModule = getTransportModule();
        m_pDynamicModule = getDynamicModule();

        print("Sequencer","Method: Transport Test\nEnd Time: %g\nOutputs: %g\n\n\n",
             configFile.getDouble("END_TIME"),
             configFile.getDouble("N_OUTPUTS"));
      printLog(std::cout);

      m_pSequencer->testTransport(*m_pDynamicModule,*m_pTransportModule,
                                  configFile.getDouble("END_TIME"),
                                  configFile.getDouble("N_OUTPUTS"));
      break;
      }
    case 2:  //Case to solve the biphasic incompressible flow
      {
        m_pTransportModule = getTransportModule();
        m_pDynamicModule = getDynamicModule();
        m_pDiff = getDiffusiveStep();
        m_pDiff->setTransport(*m_pTransportModule);


        int nDyn = configFile.getInt("N_DYN_ITERATIONS");
        int nOuts = configFile.getInt("N_OUTPUTS");
        int nOutsPerDyn = nOuts/nDyn;
        if (nOutsPerDyn ==  0)
          nOutsPerDyn = 1;
        print("Sequencer","Method: Staggered Iteration\nEnd Time: %g\nDynamic Modules Iterations: %d\nOutputs Per Dyn Iteration: %d\n",
               configFile.getDouble("END_TIME"),
               configFile.getInt("N_DYN_ITERATIONS"),
               nOutsPerDyn);
        printLog(std::cout);
        m_pSequencer->alternateIteration(*m_pDynamicModule,*m_pTransportModule,m_pDiff,
                                         configFile.getDouble("END_TIME"),
                                         configFile.getInt("N_DYN_ITERATIONS"),
                                         nOuts);
        break;
      }
    case 3:
      {
        m_pFlash = getFlash();
        m_pTransportModule = getTransportModule();
        m_pDynamicModule = getDynamicModule();
        m_pFlash->setTransport(*m_pTransportModule);
        m_pFlash->setDynamic(*m_pDynamicModule);

        m_pDiff = getDiffusiveStep();
        m_pDiff->setTransport(*m_pTransportModule);


        /* Now if this is a compressible Model the initial condition of the transport
           must be updated such that all mixture could have the volume filled buy all
           the pores. We do that buy changing the number of total moles of all the mixture
           preserving the same proportion of moles in each component previously
           defined buy the initial conditions of the transport equations 
        */
        if (typeid(*m_pFlash) != typeid(DummyFlash)) 
        {
          print("Sequencer",
              "Adjusting Initial Conditions for Compressible Model\n");
          adjustInitialConditionForCompressibleModel(*m_pDynamicModule,*dynamic_cast<ConservativeMethodForSystem*>(m_pTransportModule),*m_pFlash);
        }


        print("Sequencer","Method: Staggered Iteration with Flash calculation\n");
        print("Sequencer","End Time: %g\n"          ,configFile.getDouble("END_TIME"));
        print("Sequencer","Dynamic Iterations: %g\n",configFile.getDouble("N_DYN_ITERATIONS"));
        ConservativeMethodForSystem* pTransportSystem = dynamic_cast<ConservativeMethodForSystem*>(m_pTransportModule);


        if (configFile.isDefined("TRANSPORT_PER_DYNAMIC_STEPS"))
        {
          CompressibleDynamic *pCompDynamic;
          pCompDynamic = dynamic_cast<CompressibleDynamic*>(m_pDynamicModule);
          if (!pCompDynamic)
          {
            throw new Exception("The proportion control sequence alghorithm  demands a module that implements CompressibleDynamicBase\n");
          }
          print("Sequencer","Max Transport Steps per Cycle: %g\n",configFile.getDouble("TRANSPORT_PER_DYNAMIC_STEPS"));
          printLog(std::cout);
          m_pSequencer->alternateIterationProportionControl(*pCompDynamic,
                                                            *pTransportSystem,
                                                            *m_pFlash,m_pDiff,
                                                            configFile.getDouble("END_TIME"),
                                                            configFile.getDouble("TRANSPORT_PER_DYNAMIC_STEPS"),
                                                            configFile.getDouble("TRANSPORT_PER_DYNAMIC_STEPS_TOLERANCE"),
                                                            getDynamicTimeStep(),
                                                            configFile.getDouble("N_OUTPUTS"));

        }
        else
        {
          
          printLog(std::cout);
          m_pSequencer->alternateIteration(*m_pDynamicModule,
                                           *pTransportSystem,
                                           *m_pFlash,
                                           m_pDiff,
                                       configFile.getDouble("END_TIME"),
                                       configFile.getDouble("N_DYN_ITERATIONS"),
                                       configFile.getDouble("N_OUTPUTS"));
        }
      }
      break;
    case 4:
      {
        m_pTransportModule = getTransportModule();
        m_pDiff = getDiffusiveStep();
        m_pDiff->setTransport(*m_pTransportModule);

        int nDyn = configFile.getInt("N_DYN_ITERATIONS");
        unsigned nOuts = configFile.getInt("N_OUTPUTS");
        int nOutsPerDyn = nOuts/nDyn;
        if (nOutsPerDyn ==  0)
          nOutsPerDyn = 1;
        printLog(std::cout);
        m_pSequencer->diffusiveStepTest(*m_pDiff,                                                                        configFile.getInt("N_DYN_ITERATIONS"),
                                         configFile.getDouble("END_TIME"),
                                         nOuts);

        break;


      }
    default:
      throw new Exception("Wrong option for the SEQUENCE_ALGORITHM key\n");
      return; //Just to avoid compiling warnings
    }
    printf("\n\nSimulation Done\n");
    hdf5.close();
    return;
 
  }
  catch(Exception *e)
  {
    printf("Error");
    printBeginSection("ERROR MESSAGE");
    printf("%s\n\nQuitting.....\n",e->getError().c_str());
    abort();
  }
  
}
    

  


void MainAppMPI::printLocalTriangulation() 
{
  Point3D P0 = getLocalP0();
  Point3D P1 = getLocalP1();
  Point3D DX = getDX();
  std::vector<unsigned> el = getLocalNElems();
  print("Mesh","\
Local Triangulation:\n\
Elements: %d x %d x %d\n\
Dimensions: %g x %g %g\n\
Steps:      %g,%g,%g\n\
Domain <%g, %g, %g> - <%g,%g,%g>\n",
        el[0],el[1],el[2],
        P1[0]-P0[0],P1[1]-P0[1],P1[2]-P0[2],
        DX[0],DX[1],DX[2],
        P0[0],P0[1],P0[2],
        P1[0],P1[1],P1[2]);
  print("Mesh","Total: %u cells, %u faces\n",getMesh().numCells(),getMesh().numFaces());
}

 


 


std::string MainAppMPI::appendRankInFileName(std::string fileName) 
{

    char strHDF5[500];
    sprintf(strHDF5,"%s_%d.hdf",
            strtok( basename((char*) fileName.c_str()),"."),getRank());

    string result = dirname((char*) fileName.c_str());
    result+="/";
    result+=strHDF5;
    return result;
}


bool MainAppMPI::runningInParallel() 
{
  if (NetMPI::nProcess() == 1)
    return false;
  else
    return true;

//   int option  = configFile.getInt("DYNAMIC_MODULE");
//   int option2 = configFile.getInt("TRANSPORT_MODULE");
//   if ((option > 1 && option <=8))
//       return false;
//   else
//     return true;
}

/*
TransportBase* MainAppMPI::getTransportModule() 
{
  if (!runningInParallel())
  {
      return MainApp::getTransportModule();
  }
  this->printBeginSection("TRANSPORT MODULE");
  switch (configFile.getInt("TRANSPORT_MODULE"))
  {
  case 1:
    {
      print("Method: Lax Friedrichs with MPI\nCFL: %g\n",configFile.getDouble("MAX_CFL"));
      m_flux = getFluxFunction(); //Get the mobility function
      m_fTransBC = getTransportBC("S1"); //Get the funtion describing the boundary condition 
      m_fInit = getTransportIC("S1"); //Get the function for the initial condition
      return new LaxFriedrichsMethodMPI(*m_pMesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,*m_flux,
                                        configFile.getDouble("MAX_CFL"),
                                        configFile.getInt("MESH_OVERLAP_SIZE"));
      break;
    }
  case 2:
    throw new Exception("Method Upwind is not developed for MPI for while");
  case 3:
      print("Method: Russanov with MPI\nCFL: %g\n",configFile.getDouble("MAX_CFL"));
      m_flux = getFluxFunction(); //Get the mobility function
      m_fTransBC = getTransportBC("S1"); //Get the funtion describing the boundary condition 
      m_fInit = getTransportIC("S1"); //Get the function for the initial condition
      return new RussanovMethodMPI(*m_pMesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,*m_flux,
                                   configFile.getDouble("MAX_CFL"),
                                   configFile.getInt("MESH_OVERLAP_SIZE"));
      break;
  default:
    throw new Exception("Invalid option for TRANSPORT_MODULE selection");
    return NULL; //The return value is here just to avoid compiling warnings
  }

  
}





DynamicBase* MainAppMPI::getDynamicModule() 
{
  int option = configFile.getInt("DYNAMIC_MODULE");

  if (!runningInParallel())
  {
      return MainApp::getDynamicModule();
  }
  printBeginSection("DYNAMIC MODULE");
  switch (option)
  {
  case 1:
    {
          print("Const Velocity Module Selected:\n\t");
          return getConstVelocityDynamicModule();
   }
  case 9:
    {
      print("Hybrid finite element method with DDM\nSolver: UMFPACK\n");
      getPBoundaryCondition(m_fDP,m_fFluxBC);
      m_fK = getK();
      m_fMobT = getTotalMobility();
      m_fGravMobSrc = getGravityMobilitySource();
      int debugLevel = configFile.getInt("DYNAMIC_DEBUG_LEVEL");
      double ddmTol = configFile.getDouble("DDM_TOLERANCE");
      double b = configFile.getDouble("DDM_B_PARAMETER");
      double theta = configFile.getDouble("DDM_THETA_PARAMETER");
      int meshOverlap = configFile.getInt("MESH_OVERLAP_SIZE");


      //Permeability (Todo: Functions to get permeability)
      printf("DDM TOLERANCE: %g\nMESH_OVERLAP: %d\n\nB_PARAMETER=%g\n",ddmTol,meshOverlap,b);

      return new HybridPressionModuleUMFPACK_MPI(getMesh(),*m_fFluxBC,*m_fDP,*m_fK,*m_fMobT,*m_fGravMobSrc,
                                                 configFile.getDouble("DIM_X"),b,theta,meshOverlap,ddmTol,debugLevel,out);

    }
  case 10:
    {
      print("Hybrid finite element method with DDM and Solver: CG with AGM\n");
      getPBoundaryCondition(m_fDP,m_fFluxBC);
      m_fK = getK();
      m_fMobT = getTotalMobility();
      m_fGravMobSrc = getGravityMobilitySource();
      int debugLevel = configFile.getInt("DYNAMIC_DEBUG_LEVEL");
      double ddmTol = configFile.getDouble("DDM_TOLERANCE");
      double b = configFile.getDouble("DDM_B_PARAMETER");
      double theta = configFile.getDouble("DDM_THETA_PARAMETER");
      int meshOverlap = configFile.getInt("MESH_OVERLAP_SIZE");
      double agm_tol = configFile.getDouble("TOLERANCE");
      double agm_nIt = configFile.getInt("NUM_ITERATIONS");
      
      printf("DDM TOLERANCE: %g\nMESH_OVERLAP: %d\nAGM TOLERANCE: %g\nNUMBER OF ITERATIONS: %g\nB_PARAMETER=%g\nTHETA_PARAMETER=%g\n",
             ddmTol,meshOverlap,agm_tol,agm_nIt,b,theta);

      //Permeability (Todo: Functions to get permeability)
      return new HybridPressionModuleAGM_MPI(getMesh(),*m_fFluxBC,*m_fDP,*m_fK,*m_fMobT,*m_fGravMobSrc,
                                             getDX()[0],b,theta,meshOverlap,ddmTol,agm_tol,agm_nIt,debugLevel,out);

    }


  default:
    {
      throw new Exception("Invalid option for DYNAMIC_MODULE\n");
      return NULL;
    }

  }
}
*/

---