AdlibCoupledSolver3.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ludovic Noels, Ralf Deiterding

#ifndef ADLIB_COUPLED_SOLVER3_H
#define ADLIB_COUPLED_SOLVER3_H

#define FACTOR_ADLIB 1.0e2
#include <cfloat>
#include <cmath>

#include "Solver.h"

// adlib header files
#include "mesher3d/mesher3d.h"
#include "mechanics/mechanics.h"
#include "materials/matlib.h"
#include "mechanics/readmesh.h"
#include "mesher3d/importmesh.h"
#include "mechanics/initialconditions.h"
#include "external_coupling_3d/external_coupling_3d.h"
#include "external_coupling_3d/boundaryconditions.h"
#include "mpi_adlib/mpi_adlib.h"
#include "external_coupling/AdlibTiming.h"

// There is in fact nothing dimension-specific in this class. Specializing with dim=3 is only
// necessary, because the directories adlib/external_coupling_2d and adlib/mesher_2d are in a 
// very old state and need to be brought up to date. 

template <>
class AdlibCoupledSolver<3> : public Solver {
  typedef Solver base;
public:
  AdlibCoupledSolver(MPI_Comm solidComm) : base(), SubIterations(1), _olog(0) {

    setCurrentTime(0.);
    setCurrentTimeStep(0.);
    setCurrentStepCount(0);
    renameOldCheckpointedFiles=false;

    // initialize solid processor logfile
    MPI_Comm_rank(solidComm, &myRank);
#ifdef DEBUG_PRINT
    std::ostringstream obuf;
    obuf << "S" << myRank << ".log" << static_cast<char>(0);
    oflog.open(obuf.str().c_str());     
    _olog = new std::ostream(oflog.rdbuf());
#endif

    // initialize adlib solid mechanics communicator
    adlib::AdlibSetMechanicsCommunicator(solidComm);

    fsBoundary=NULL;

    scaleBoundary = false;
    xCenter = 0.; 
    yCenter = 0.; 
    zCenter = 0.;
    ratioXPos = 1.; 
    ratioXNeg = 1.;
    ratioYPos = 1.; 
    ratioYNeg = 1.;
    ratioZPos = 1.; 
    ratioZNeg = 1.;

  }

  ~AdlibCoupledSolver() {       
    if (_olog) delete _olog;
    if (fsBoundary) delete fsBoundary;
  }

  //******************************************************************************
  // Abstract class interface
  //******************************************************************************
  virtual void sendBoundaryReceivePressure() = 0;  
  //******************************************************************************

  // Read input files, etc.
  virtual bool setup() {

#ifdef DEBUG_PRINT
    ( log() << "Enter setup\n" ).flush();
#endif
#ifdef DEBUG
    // set diagnostic output level
    adlib::prt = 1;
#else
    adlib::prt = 0;
#endif

    // read adlib dimensions
    adlib::dimensions("dimensions.dat");
 
    // read input control parameters from inp
    std::ifstream in("inp");
    in.getline(meshFile,80);
    in.getline(listMatFile,80);
    in.getline(listTypeFile,80);
    in.getline(listBoundBBFile,80);
    in.getline(listScaleBoundaryFile,80);
    in.getline(materialFile,80);
    in >> time_factor;
    in >> subdivisions;
    in >> SubIterations;
    in >> gamma_newmark;
    in.getline(outdir,80);
    in.close();

    if (gamma_newmark < 0.5)
      gamma_newmark = 0.5;
    if (gamma_newmark > 1.5)
      gamma_newmark = 1.5;
    if (SubIterations<=0) 
      SubIterations=1;

#ifdef DEBUG_PRINT
    ( log() << "Gamma Newmark = " << gamma_newmark << "\n").flush();
    ( log() << "Spectral radius at bifurcation = " << 1.-4.*(gamma_newmark-0.5)/(gamma_newmark+0.5)/(gamma_newmark+0.5) << "\n").flush();
#endif

    ADLIB_START_WATCH_WHOLE

    return true;
  }

  virtual void finish() {
    ADLIB_END_WATCH_WHOLE
    AdlibTiming::collect(*adlib::AdlibGetMechanicsCommunicator());
#ifdef DEBUG_PRINT
    AdlibTiming::print_local(log());
#endif
    if (adlib::AdlibMechanicsProcId() == 0)
      AdlibTiming::print(std::cout);
  }

  virtual bool readMeshFiles(bool isARestart)
  {
    if(!isARestart) {
      // read in solid mesh and material data
      double * coordinates=NULL;
      int * connectivities=NULL;
      int * elementMaterials=NULL;
      int nbNodesElementRead=0;
        
      adlib::readMeshFromTextFile(meshFile,
                                  &nbNodesElementRead,
                                  &adlib::nodes,
                                  &adlib::elements,
                                  &coordinates,
                                  &connectivities);

      if (strcmp(listMatFile,"default")) 
        adlib::readElementMaterialsFromTextFile(listMatFile,
                                                adlib::elements,
                                                &elementMaterials);
        
      if (strcmp(listTypeFile,"default"))
        adlib::readElementTypesFromTextFile(listTypeFile,
                                            adlib::elements,
                                            &adlib::element_type);
        
      adlib::importMesh(nbNodesElementRead,
                        adlib::nodes,
                        adlib::elements,
                        coordinates,
                        connectivities,
                        elementMaterials,
                        0,
                        NULL,
                        NULL);

      // initialize boundary and forces to zero by hand
      std::memset(adlib::boundary, 0, sizeof(int)*(adlib::nodes));
      std::memset(adlib::forces, 0, 
                  sizeof(double)*(adlib::nodes)*(adlib::dof_node));

      // output mesh statistics and dump imported mesh
      //adlib::TetrahedraStatistics();
      //adlib::tecplot("both","mesh","mesh",0,0,"initialMesh",1.,"ascii");

#ifdef DEBUG_PRINT
      ( log() << "Total number of elements:" << adlib::elements << "\n").flush();
#endif

      // partition imported solid mesh
      if (nProcessors > 1)
        adlib::PartitionMesh(nProcessors,myPId,subdivisions,"metis");

#ifdef DEBUG_PRINT
      ( log() << "Partioning in " << adlib::elements << "\n" ).flush();
#endif
    }

    // adlib initialization tasks
    adlib::profile();
    adlib::shape_functions();

    if(!isARestart) 
      adlib::initialize_strains();

    // read material model parameters and assemble mass
    adlib::ReadMaterialClasses(materialFile);
    adlib::assemble("mass","deformation");
    
    // read boundary conditions
    if (strcmp(listBoundBBFile,"default"))
      adlib::readAndImportFixationInBoxFromTextFile(listBoundBBFile);

    // read scaling boundary
    if (strcmp(listScaleBoundaryFile,"default"))
      readScalingBoundaryFile();

    // boundary
    if(!isARestart) {
      // generate solid/fluid boundary
      if(!fsBoundary) fsBoundary = new adlib::Boundary;
      *fsBoundary = adlib::AdlibBoundaryGenerate(adlib::solidFluid);

#ifdef DEBUG_PRINT
      ( log() << "Solid fluid boundary created\n" ).flush();
#endif

      // dump solid/fluid boundary
      char name[80];
      std::sprintf(name,"fs-boundary-%03d.tec",myRank);
      adlib::AdlibBoundaryTecplot(name, fsBoundary);
    }
    return true;
  }

  virtual void Initialize(double& t, double& dt) {
    ADLIB_START_WATCH
    readMeshFiles(false);
    sendBoundaryReceivePressure();
    dt = NextTimeStep()*SubIterations;
    t = CurrentTime();
    ADLIB_END_WATCH(INITIALIZATION)
  }

  virtual void Advance(double& t, double& dt) {

    double tf = t+dt;
    dt /= SubIterations;

    sendBoundaryReceivePressure();

#ifdef DEBUG_PRINT
    ( log() << "--- Iteration: " << CurrentStepCount() 
            << " t = " << t 
            << " dt = " << dt <<"\n" ).flush();
#endif 

    while (tf-t > DBL_EPSILON*FACTOR_ADLIB) {   
      Update(dt);
      t = CurrentTime();
      dt = NextTimeStep();
      if (dt > tf-t) dt = tf-t;
    }

    setCurrentStepCount(CurrentStepCount() + SubIterations);
    dt = NextTimeStep()*SubIterations;
    t = tf;
    setCurrentTime(tf);
  }

  // Numerical Update
  virtual void Update(const double& dt) {

    setCurrentTimeStep(dt);
    ADLIB_START_WATCH
      //adlib::explicitIntegration("cdif","pred");
      adlib::newmark("expl",0.,gamma_newmark);
    ADLIB_END_WATCH(INTEGRATION)
    ADLIB_START_WATCH
      adlib::assemble("explicit","deformation");
    ADLIB_END_WATCH(ASSEMBLE)
    ADLIB_START_WATCH
      //adlib::explicitIntegration("cdif","corr");
      adlib::newmark("accel",0.,gamma_newmark);
    ADLIB_END_WATCH(INTEGRATION)
    setCurrentTime(CurrentTime()+dt);
  }

  const double& CurrentTime() const { return adlib::time_n; }
  void setCurrentTime(double _time) { adlib::time_n = _time; }

  const double& CurrentTimeStep() const { return adlib::time_step; }
  void setCurrentTimeStep(double _dt) { adlib::time_step = _dt; }

  const int& CurrentStepCount() const { return adlib::step_number; }
  void setCurrentStepCount(int _step) { adlib::step_number = _step; }

  double NextTimeStep() const {
      return time_factor*(1./(0.5+gamma_newmark))*adlib::StableTimeStep("deformation");
  }
  
  virtual void Output() {
    ADLIB_START_WATCH

    // Write data files
    int sN=CurrentStepCount()/SubIterations;
    // boundary

    char name[80];
    std::sprintf(name,"solidBoundary-%07d-%03d.tec",sN,myRank);
    adlib::AdlibBoundaryTecplot(name, fsBoundary);

    // internal
    std::sprintf(name,"internal-%07d-%03d",sN,myRank); 
    char timeChar[16];
    std::sprintf(timeChar,"%f*s",adlib::time_n);
    adlib::tecplot("both",name,"internal",0,1,timeChar,1.,"ascii");

    // stresses
    std::sprintf(name,"stresses-%07d-%03d",sN,myRank);
    adlib::tecplot("both",name,"stresses",0,1,timeChar,1.,"ascii");

    if (sN == 0)
      {
        std::sprintf(name,"material-%03d",myRank);          
        adlib::materialOutput(name);
      }

    ADLIB_END_WATCH(OUTPUT)
  }

  virtual int NSteps() {
    return CurrentStepCount()/SubIterations;
  }

  virtual void Restart(double& t, double& dt) {

#ifdef DEBUG_PRINT  
    ( log() << "Entering restart. \n" ).flush();
#endif 

    // Read solid mesh 
    char nameRestartMesh[80];
    sprintf(nameRestartMesh,"adlib-mesh.res.%03d.cp",myPId);
    adlib::RestartFileRead(nameRestartMesh,1); 
  
    // Read parallel solid mechanics communicator data 
    if (nProcessors > 1) {
      char nameRestartCom[80];
      sprintf(nameRestartCom,"adlib-com.res.%03d.cp",myPId); 
      adlib::PMechRestartRead(myPId,nameRestartCom); 
    }

    // read fluid solid boundary
    char nameRestartBound[80];
    sprintf(nameRestartBound,"adlib-bound.res.%03d.cp",myPId);
    if(!fsBoundary) fsBoundary = new adlib::Boundary;
    *fsBoundary=adlib::AdlibBoundaryRead(nameRestartBound);

    readMeshFiles(true);

    //proceed
    t = CurrentTime();
    dt = NextTimeStep()*SubIterations;
    sendBoundaryReceivePressure();

    // boundary
    char name[80];
    std::sprintf(name,"solidBoundaryRestart-%03d.tec",myRank);
    adlib::AdlibBoundaryTecplot(name, fsBoundary);

    // internal
    std::sprintf(name,"internalRestart-%03d",myRank); 
    char timeChar[16];
    std::sprintf(timeChar,"%f*s",adlib::time_n);
    adlib::tecplot("both",name,"internal",0,1,timeChar,1.,"ascii");

  }

  virtual void Checkpointing() {
    ADLIB_START_WATCH

    // Write solid mesh
    char nameRestartMesh[80];
    sprintf(nameRestartMesh,"adlib-mesh.res.%03d.cp",myPId);
    if(renameOldCheckpointedFiles)
    {
      char oldNameRestartMesh[80];
      sprintf(oldNameRestartMesh,".adlib-mesh.res.%03d.cp",myPId);
      rename(nameRestartMesh,oldNameRestartMesh);
    }

    adlib::RestartFileWrite(nameRestartMesh,1); //or 0?

    // Write parallel solid mechanics communicator data
    if (nProcessors > 1) {
      char nameRestartCom[80];
      sprintf(nameRestartCom,"adlib-com.res.%03d.cp",myPId);
      if(renameOldCheckpointedFiles)
      {
        char oldNameRestartCom[80];
        sprintf(oldNameRestartCom,".adlib-com.res.%03d.cp",myPId);
        rename(nameRestartCom,oldNameRestartCom);
      }
      adlib::PMechRestartWrite(nameRestartCom);
    }
    // Write fluid solid boundary
    char nameRestartBound[80];
    sprintf(nameRestartBound,"adlib-bound.res.%03d.cp",myPId);
    if(renameOldCheckpointedFiles)
    {
      char oldNameRestartBound[80];
      sprintf(oldNameRestartBound,".adlib-bound.res.%03d.cp",myPId);
      rename(nameRestartBound,oldNameRestartBound);
    }
    adlib::AdlibBoundaryWrite(nameRestartBound, fsBoundary);

    renameOldCheckpointedFiles=true;

    ADLIB_END_WATCH(OUTPUT)
  }

  virtual void readScalingBoundaryFile()
  {
    scaleBoundary = true;
    std::ifstream in(listScaleBoundaryFile);
    in >> xCenter;
    in >> ratioXPos;
    in >> ratioXNeg;
    in >> yCenter;
    in >> ratioYPos;
    in >> ratioYNeg;
    in >> zCenter;
    in >> ratioZPos;
    in >> ratioZNeg;
    
#ifdef DEBUG_PRINT
      ( log() << "Scaling boundary:\n" ).flush();
      ( log() << "  around x = " << xCenter << " by a factor " << ratioXPos 
        << " for x > " << xCenter << " and by a factor " << ratioXNeg 
        << " for x < " << xCenter << "\n" ).flush();
      ( log() << "  around y = " << yCenter << " by a factor " << ratioYPos 
        << " for y > " << yCenter << " and by a factor " << ratioYNeg 
        << " for y < " << yCenter << "\n" ).flush();
      ( log() << "  around z = " << zCenter << " by a factor " << ratioZPos 
        << " for z > " << zCenter << " and by a factor " << ratioZNeg 
        << " for z < " << zCenter << "\n" ).flush();
#endif

  }


  inline const adlib::Boundary& Boundary() const { return *fsBoundary; }
  inline adlib::Boundary& Boundary() { return *fsBoundary; }

  inline std::ostream& log() { return *_olog; }

protected:
  int SubIterations;
  std::ostream* _olog;
  std::ofstream oflog; 
  int myRank;

  // input parameters
  char meshFile[80]; // solid mesh descriptor file name
  char listMatFile[80]; // solid mesh element materials file name
  char listTypeFile[80]; // solid mesh element type file name
  char listBoundBBFile[80]; // boundary conditions file name (bounding box)
  char listScaleBoundaryFile[80]; // scaling boundary file name
  char materialFile[80]; // solid material model parameter file name
  int subdivisions; // levels of subdivision of initial solid mesh
  char outdir[80]; // name of output directory
  bool renameOldCheckpointedFiles;

  // scaling boundary parameters
  bool scaleBoundary;
  double xCenter, yCenter, zCenter;
  double ratioXPos, ratioXNeg;
  double ratioYPos, ratioYNeg;
  double ratioZPos, ratioZNeg;


private:

  double time_factor;
  double gamma_newmark; 
  adlib::Boundary* fsBoundary;
};

#endif

---