amroc/amr/AMRELCGFMSolver.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@cacr.caltech.edu

#ifndef AMROC_ELC_GFMSOLVER_H
#define AMROC_ELC_GFMSOLVER_H

#include "AMRCoupledGFMSolver.h"
#include "AMRGFMInterpolation.h"
#include "CPTLevelSet.h"
#include "Timing.h"
#include "elc/EulerianCommBoundary.h"
#include "elc/EulerianCommShell.h"
#include <vector>


template <class VectorType, class FixupType, class FlagType, int dim>
class AMRELCGFMSolver : public AMRCoupledGFMSolver<VectorType,FixupType,FlagType,dim> {
  typedef AMRCoupledGFMSolver<VectorType,FixupType,FlagType,dim> base;
  typedef typename VectorType::InternalDataType DataType;
public:
  typedef typename base::integrator_type integrator_type;
  typedef typename base::initial_condition_type initial_condition_type;
  typedef typename base::boundary_conditions_type boundary_conditions_type;
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  

  typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,dim> interpolation_type;
  typedef typename interpolation_type::point_type point_type;
  typedef typename elc::EulerianComm<dim,DataType> eul_comm_type;
  typedef typename elc::EulerianCommShell<dim,DataType>    eul_comm_pressure_on_face_type;
  typedef typename elc::EulerianCommBoundary<dim,DataType> eul_comm_pressure_on_node_type;
  typedef typename geom::BBox<dim,DataType> bbox_type;
  typedef CPTLevelSet<DataType,dim> cpt_type;

  AMRELCGFMSolver(integrator_type& integ, 
                  initial_condition_type& init,
                  boundary_conditions_type& bc) : base(integ, init, bc), _CoupleGFM(-1),
    _SendDataLocation(0), _ThinStructure(0), _LengthConversionFactor(1.), _PressureConversionFactor(1.),
    _SafeCFL(1.) {      
    _Interpolation = new interpolation_type(*this);
    _eulComm = (eul_comm_type*) 0; 
  }

  virtual ~AMRELCGFMSolver() {
    delete _Interpolation;
    if (_eulComm) delete _eulComm;
  }
    
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
    RegisterAt(base::LocCtrl,"DataLocation",_SendDataLocation);
    RegisterAt(base::LocCtrl,"ThinStructure",_ThinStructure);
    RegisterAt(base::LocCtrl,"LengthFactor",_LengthConversionFactor);
    RegisterAt(base::LocCtrl,"PressureFactor",_PressureConversionFactor);
    RegisterAt(base::LocCtrl,"SafeCFL",_SafeCFL);
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    base::register_at(Ctrl);
  }

  virtual void finish() {
    if (_eulComm) {
      delete _eulComm;
      _eulComm = (eul_comm_type*) 0;
    }
    base::finish();
  }

  virtual void SetupData() {
    base::SetupData();
    _Interpolation->SetupData(base::PGH(), base::NGhosts());
    if (_SendDataLocation!=0) _SendDataLocation=1;
  }

  void SetupInterComm(const int solid_nodes, const int solid_root) {
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    if (SendDataLocation()==0)
      _eulComm = new eul_comm_pressure_on_face_type(MPI_COMM_WORLD, base::Comm(), solid_nodes,
                                                    solid_root, elc::GlobalIdentifiers);
    else
      _eulComm = new eul_comm_pressure_on_node_type(MPI_COMM_WORLD, base::Comm(), solid_nodes,
                                                    solid_root, elc::GlobalIdentifiers);
#ifdef DEBUG_PRINT_ELC
    ( comm_service::log() << "*** EulerianComm created: SolidNodes=" << solid_nodes 
      << " SolidRoot=" << solid_root << "\n" ).flush();
#endif
  }

  double CFLGather(double cfl) {
#ifdef DAGH_NO_MPI
#else
    if (comm_service::dce() && comm_service::proc_num() > 1) { 
      
      int num = comm_service::proc_num();
      
      int R;
      int size = sizeof(double);
      void *values = (void *) new char[size*num];
      
      START_WATCH;
      R = MPI_Allgather(&cfl, size, MPI_BYTE, values, size, MPI_BYTE, 
                        comm_service::comm(base::U().GF_Id()));
      END_WATCH(GATHERCFL);
      if ( MPI_SUCCESS != R ) 
        comm_service::error_die( "AMRELCGFMSolver::cfl", "MPI_Allgather", R );
      
      for (int i=0;i<num;i++) {
        double tmax = *((double *)values + i);
        cfl = cfl > tmax ? cfl : tmax;
      }
      if (values) delete [] (char *)values;
    }
#endif
    return cfl;
  }

  // Temporary solution until local time step adjustment is implemented correctly
  // in a future version. 
  virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level, 
                                double t, double dt, bool DoFixup, double tc, 
                                const int which) {
    double cfl_orig = base::IntegrateLevel(u,Time,Level,t,dt,DoFixup,tc,which);
    if (_SafeCFL > 0.) {
      double cfl_max = CFLGather(cfl_orig);
      int NSteps = int(cfl_max / _SafeCFL) + 1;
      while (cfl_max > _SafeCFL) {
        double tl = t;
        double dtl = dt/NSteps;
#ifdef DEBUG_PRINT_AMRSOLVER
        ( comm_service::log() << "IntegrateLevel(): Level " << Level 
          << "  to be repeated with " << NSteps << "  steps. Global CFL = " << cfl_max << "\n" ).flush();
#endif
        cfl_max = 0.;
        for (int i=1; i<=NSteps; i++) {
          if (i==NSteps) dtl=t+dt-tl;
          double cfl = base::IntegrateLevel(u,Time,Level,tl,dtl,DoFixup,tc,which);
          cfl_max = cfl_max > cfl ? cfl_max : cfl;
          tl += dtl;
        }   
        cfl_max = CFLGather(cfl_max);
#ifdef DEBUG_PRINT_AMRSOLVER
        ( comm_service::log() << "IntegrateLevel(): Level " << Level 
          << "  repeated with " << NSteps << "  steps. New global CFL = " << cfl_max << "\n" ).flush();
#endif
        NSteps++;
      }
    }
    return cfl_orig;
  }

  virtual void SendBoundaryData() {
#ifdef DEBUG_PRINT
      ( comm_service::log() << "*** SendBoundaryData()\n" ).flush();
#endif
      int Time = CurrentTime(base::GH(),base::CouplingLevel);

      // Linear interpolation of coarser level data that is not at the same time level
      for (register int l=0; l<=FineLevel(base::GH()); l++) {
        int TimeC = CurrentTime(base::GH(),l);
        int TStepC = TimeStep(base::Work(),l);
        if (TimeC < Time) {
          double frac = 1.0*(Time-TimeC)/TStepC;
          forall (base::Work(),TimeC,l,c)   
            base::Work()(TimeC,l,c).lin_interp(base::Work()(TimeC,l,c),(1.0-frac),
                                               base::Work()(TimeC+TStepC,l,c),frac);
          end_forall
        }
        else if (TimeC+TStepC == Time) {
          forall (base::Work(),TimeC,l,c)   
            base::Work()(TimeC,l,c).copy(base::Work()(TimeC+TStepC,l,c));
          end_forall
        }
      }
    START_WATCH
#ifdef DEBUG_PRINT
      ( comm_service::log() << "*** Geometry construction()\n" ).flush();
#endif
      DataType distance = -base::GFM(CoupleGFM()).Boundary().Cutoff();
    
      DataType dxmin = 1.e37;
      forall(base::U(),Time,base::CouplingLevel,c)
        DCoords dx = base::GH().worldStep(base::U()(Time,base::CouplingLevel,c).stepsize());
        for (register int d=0; d<base::Dim(); d++) 
          if (dx(d)<dxmin) dxmin=dx(d);
      end_forall 
      if (distance>dxmin*base::NGhosts()) {
        int ghn = static_cast<int>(distance/dxmin + 1);
        std::cerr << "\n\nToo few ghostcells for Cutoff=" << -distance 
                  << ". Set GhostCells>=" << ghn << " or set Cutoff<=" 
                  << dxmin*base::NGhosts() << ".\n";
        std::cerr << "Aborting programm...\n" << std::flush;
        std::exit(-1);
      }

      _eulComm->initializePressure();

      register int n;
      int Npoints;

      double *press_data = _eulComm->getPressuresData();
      const int* con = _eulComm->getConnectivitiesData();
      const point_type *pos = reinterpret_cast<const point_type *>(_eulComm->getPositionsData());
      if (_ThinStructure>0 || distance>1.e-12)
        _eulComm->computeFaceNormals();
      if (_ThinStructure>0 || SendDataLocation()!=1 || distance>1.e-12)
        _eulComm->computeFaceCentroids();
      const point_type *norm = reinterpret_cast<const point_type *>(_eulComm->getFaceNormalsData());
      const point_type *cent = reinterpret_cast<const point_type *>(_eulComm->getFaceCentroidsData());

      if (SendDataLocation()==1)
        Npoints = _eulComm->getNumberOfNodes();
      else
        Npoints = _eulComm->getNumberOfFaces();

      if (_ThinStructure>0) {
        point_type* positions = new point_type[2*Npoints];
        point_type* centroids = new point_type[2*Npoints];
        if (SendDataLocation()==1) {
          DataType* face_count = new DataType[Npoints];
          point_type* point_normal = new point_type[Npoints];
          START_WATCH
            for (n=0; n<Npoints; n++) {
              face_count[n] = static_cast<DataType>(0.);
              point_normal[n] = static_cast<point_type>(0.);
            }
              
            int Nfaces = _eulComm->getNumberOfFaces();
            for (n=0; n<Nfaces; n++) 
              for (register int d=0; d<dim; d++) {
                point_normal[con[dim*n+d]] += norm[n]; 
                face_count[con[dim*n+d]] += static_cast<DataType>(1.);
              }

            for (n=0; n<Npoints; n++) {
              point_normal[n] /= face_count[n];
              positions[2*n]   = pos[n]-distance*point_normal[n];
              positions[2*n+1] = pos[n]+distance*point_normal[n];
              centroids[2*n+1] = centroids[2*n] = pos[n];
            }

          END_WATCH(FLUID_CPL_ELC_GEOMETRY)
          delete [] face_count;
          delete [] point_normal;
        }
        else {
          START_WATCH
            for (n=0; n<Npoints; n++) {
              positions[2*n]   = cent[n]-distance*norm[n];
              positions[2*n+1] = cent[n]+distance*norm[n];
              centroids[2*n+1] = centroids[2*n] = cent[n];
            }
          END_WATCH(FLUID_CPL_ELC_GEOMETRY)
        }
#ifdef DEBUG_PRINT
        ( comm_service::log() << "*** Pressure interpolation()\n" ).flush();
#endif
        DataType* press_int = new DataType[2*Npoints];

        START_WATCH
          int TimeZero = CurrentTime(base::GH(),0);
          int Lev=0, Np=2*Npoints;
          bool* lpos = new bool[Np];
          bool* lcent = new bool[Np];
          _Interpolation->LocalPoints(base::Work(),TimeZero,Lev,Np,positions,lpos);
          _Interpolation->LocalPoints(base::Work(),TimeZero,Lev,Np,centroids,lcent);
          int nlpos=0, nppos=0; 
          for (n=0; n<Np; n++)
            if (lcent[n]) {
              if (lpos[n]) nlpos++;
              else nppos++;
            }
          int lm, pm;
          point_type* lpositions = new point_type[nlpos];
          point_type* ppositions = new point_type[nppos];
          DataType* lpress_int = new DataType[nlpos];
          DataType* ppress_int = new DataType[nppos];
          for (n=0, lm=0, pm=0; n<Np; n++)
            if (lcent[n]) {
              if (lpos[n]) lpositions[lm++] = positions[n];
              else ppositions[pm++] = positions[n];
            }
          _Interpolation->PointsValues(base::Work(),base::t[0],nlpos,lpositions,lpress_int);
          _Interpolation->PointsValuesPar(base::Work(),base::t[0],nppos,ppositions,ppress_int);
          for (n=0, lm=0, pm=0; n<Np; n++)
            if (lcent[n]) 
              if (lpos[n]) press_int[n] = lpress_int[lm++];
              else press_int[n] = ppress_int[pm++];
            else press_int[n] = _Interpolation->ErrorValue();
          delete [] lpositions;
          delete [] ppositions;
          delete [] lpress_int;
          delete [] ppress_int;
          delete [] lpos;
          delete [] lcent;
        END_WATCH(FLUID_CPL_INTERPOLATE)
          
        for (n=0; n<Npoints; n++)
          if (press_int[2*n]==_Interpolation->ErrorValue() || 
              press_int[2*n+1]==_Interpolation->ErrorValue()) {
            press_data[n] = std::numeric_limits<DataType>::max();
#ifdef DEBUG_PRINT
            if (press_int[2*n]!=press_int[2*n+1] && base::GFM(CoupleGFM()).Boundary().GetVerbose())
              ( comm_service::log() << "*** Thin structure interpolation error for center at " 
                << centroids[2*n] << "   " << positions[2*n] << ": " << press_int[2*n] << "   " 
                << positions[2*n+1] << ": " << press_int[2*n+1] << "\n" ).flush();
#endif
          }
          else
            press_data[n] = press_int[2*n]-press_int[2*n+1];
        
        delete [] positions;
        delete [] centroids;
        delete [] press_int;
      }
      else {
        if (SendDataLocation()==1) {
          if (distance>1.e-12) {
            int Nfaces = _eulComm->getNumberOfFaces();
            DataType* face_count = new DataType[Npoints];
            point_type* point_normal = new point_type[Npoints];
            START_WATCH
              for (n=0; n<Npoints; n++) {
                face_count[n] = static_cast<DataType>(0.);
                point_normal[n] = static_cast<point_type>(0.);
              }
              
              for (n=0; n<Nfaces; n++) 
                for (register int d=0; d<dim; d++) {
                  point_normal[con[dim*n+d]] += norm[n]; 
                  face_count[con[dim*n+d]] += static_cast<DataType>(1.);
                }

              point_type* positions = new point_type[Npoints];
              for (n=0; n<Npoints; n++) {
                point_normal[n] /= face_count[n];
                positions[n] = pos[n]+distance*point_normal[n];
              }
            END_WATCH(FLUID_CPL_ELC_GEOMETRY)
            delete [] face_count;
            delete [] point_normal;

            START_WATCH
              _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,positions,press_data);
            END_WATCH(FLUID_CPL_INTERPOLATE)

            delete [] positions;
          }
          else {
            START_WATCH
              _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,pos,press_data);
            END_WATCH(FLUID_CPL_INTERPOLATE)
          }
        }
        else {
          point_type* positions = new point_type[Npoints];
          START_WATCH
            for (n=0; n<Npoints; n++) { 
              positions[n] = cent[n];
              if (distance>1.e-12) 
                positions[n] += distance*norm[n];
            }   
          END_WATCH(FLUID_CPL_ELC_GEOMETRY)
        
          START_WATCH
            _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,positions,press_data);
          END_WATCH(FLUID_CPL_INTERPOLATE)

          delete [] positions;
        }
        for (n=0; n<Npoints; n++) 
          if (press_data[n]==_Interpolation->ErrorValue()) 
            press_data[n] = std::numeric_limits<DataType>::max();
      }
      
#ifdef DEBUG_PRINT_ELC
      ( comm_service::log() << "*** ELC::SendPressure" ).flush();
#endif
    
      if (PressureConversionFactor()!=1.) {
        for (n=0; n<Npoints; n++) 
          if (press_data[n]!=std::numeric_limits<DataType>::max())
            press_data[n] *= static_cast<DataType>(PressureConversionFactor());
      }

      ModifyELCSendData(_eulComm);

      START_WATCH
        _eulComm->sendPressure();
        _eulComm->waitForPressure();
      END_WATCH(ELC_SENDPRESSURE)

#ifdef DEBUG_PRINT_ELC
      ( comm_service::log() << " done.\n" ).flush();
#endif
#ifdef DEBUG_PRINT_ELC_VALUES
      ( comm_service::log() << "*** Values passed to ELC ***\n" ).flush();
      for (n=0; n<Npoints; n++) 
        ( comm_service::log() << press_data[n] << "\n" ).flush();    
#endif
    END_WATCH(FLUID_CPL_SEND_OVERHEAD)
  }

  virtual void PostReceiveBoundaryData(bool FullDomain=false) {
    START_WATCH
      BBox wholebb;

      if (!FullDomain) {
        int Level = 0;
        int Time = CurrentTime(base::GH(),Level);
        forall (base::U(),Time,Level,c)      
          wholebb += base::U().boundingbbox(Time,Level,c);
        end_forall
      }
      else 
        wholebb = base::GH().wholebbox();

      DCoords lbcorner = base::GH().worldCoords(wholebb.lower(), wholebb.stepsize());
      DCoords ubcorner = base::GH().worldCoords(wholebb.upper()+wholebb.stepsize(), 
                                                wholebb.stepsize());
      if (LengthConversionFactor()!=1.) {
        lbcorner /= LengthConversionFactor();
        ubcorner /= LengthConversionFactor();
      }
      bbox_type lregion(lbcorner(), ubcorner());
            
#ifdef DEBUG_PRINT_ELC
      ( comm_service::log() << "*** ELC::PostReceiveBoundaryData(), " << lregion ).flush();
#endif

      START_WATCH
        _eulComm->receiveMesh(lregion);
      END_WATCH(ELC_RECEIVEBOUNDARY)

#ifdef DEBUG_PRINT_ELC
      ( comm_service::log() << " done.\n" ).flush();
#endif
    END_WATCH(FLUID_CPL_RECEIVE_OVERHEAD)
  }

  virtual void WaitReceiveBoundaryData() {
    START_WATCH
#ifdef DEBUG_PRINT_ELC
      ( comm_service::log() << "*** ELC::WaitReceiveBoundaryData()" ).flush();
#endif

      START_WATCH
        _eulComm->waitForMesh();
      END_WATCH(ELC_RECEIVEBOUNDARY)

      ModifyELCReceiveData(_eulComm);

#ifdef DEBUG_PRINT_ELC
      ( comm_service::log() << " done.\n" ).flush();
#endif
      int Npoints = _eulComm->getNumberOfNodes();
      point_type* positions = (point_type *)(_eulComm->getPositionsData());
      point_type* velocities = (point_type *)(_eulComm->getVelocitiesData());
      register int n;
      if (LengthConversionFactor()!=1.) {
        for (n=0; n<Npoints; n++) {
          positions[n] *= static_cast<DataType>(LengthConversionFactor());
          velocities[n] *= static_cast<DataType>(LengthConversionFactor());
        }
      }
#ifdef DEBUG_PRINT_ELC_VALUES
      ( comm_service::log() << "*** Positions received from ELC ***\n" ).flush();
      for (n=0; n<Npoints; n++) 
        ( comm_service::log() << positions[n] << "\n" ).flush();    
#endif
    END_WATCH(FLUID_CPL_RECEIVE_OVERHEAD)

      (reinterpret_cast<cpt_type* >(base::GFM(CoupleGFM()).GetLevelSetP()))->SetBrep(_eulComm->getNumberOfNodes(),
                 _eulComm->getPositionsData(), _eulComm->getNumberOfFaces(),
                 _eulComm->getConnectivitiesData());
  }

  void ExtractBoundaryVelocities(const BBox& bb, const int& nc, 
                                 const point_type* xc, point_type* retvel) {
    START_WATCH
      int nodes = _eulComm->getNumberOfNodes();
      const point_type *positions = reinterpret_cast<const point_type *>(_eulComm->getPositionsData());
      const point_type *velocities = reinterpret_cast<const point_type *>(_eulComm->getVelocitiesData());

#ifdef DEBUG_PRINT_ELC_VALUES
      ( comm_service::log() << "*** Velocities constructed in " << bb << " ***\n" ).flush();
#endif

      // We are looking here only in a Box, which is not necessarily enough
      // temporary solution until CPT is used for that
      DCoords lbcorner = base::GH().worldCoords(bb.lower(),bb.stepsize());
      DCoords ubcorner = base::GH().worldCoords(bb.upper()+bb.stepsize(),bb.stepsize());
      point_type lb, ub;
      register int m, n;
      for (m=0; m<base::Dim(); m++) {
        lb(m) = lbcorner(m); ub(m) = ubcorner(m);
      }

      std::vector<int> pos_list;
      for (m=0; m<nodes; m++)
        if (!(lb>positions[m] || ub<positions[m]))
          pos_list.push_back(m);    
      int s = pos_list.size();

      if (s>0) {
        double min_dist, dist;
        register int min_point;
        for (n=0; n<nc; n++) {
          min_dist = Abs(xc[n]-positions[pos_list[0]]);
          min_point = pos_list[0];
          for (m=1; m<s; m++) {
            dist = Abs(xc[n]-positions[pos_list[m]]);
            if (dist<min_dist) {
              min_dist = dist;
              min_point = pos_list[m];
            }
          }
          retvel[n] = velocities[min_point];      
        }
      }
      else 
        for (n=0; n<nc; n++) 
          retvel[n] = 0.;
   
#ifdef DEBUG_PRINT_ELC_VALUES
      for (n=0; n<nc; n++) 
        ( comm_service::log() << xc[n] << "   " << retvel[n] << "\n" ).flush();
#endif
    END_WATCH(FLUID_CPL_VELOCITY_SEARCH)
  }

  virtual void ModifyELCReceiveData(eul_comm_type* eulComm) {}
  virtual void ModifyELCSendData(eul_comm_type* eulComm) {}

  inline void SetCoupleGFM(int gfm) { _CoupleGFM = gfm; }
  inline const int& CoupleGFM() const { return _CoupleGFM; }
  inline void SetSendDataLocation(int bd) { SendDataLocation = bd; }
  inline const int& SendDataLocation() const { return _SendDataLocation; }
  inline void SetLengthConversionFactor(DataType l) { _LengthConversionFactor = l; }
  inline const DataType& LengthConversionFactor() const { return _LengthConversionFactor; }
  inline void SetPressureConversionFactor(DataType p) { _PressureConversionFactor = p; }
  inline const DataType& PressureConversionFactor() const { return _PressureConversionFactor; }

protected:
  interpolation_type* _Interpolation;
  eul_comm_type* _eulComm;
  int _CoupleGFM;
  int rank;
  int _SendDataLocation, _ThinStructure;
  DataType _LengthConversionFactor, _PressureConversionFactor;
  double _SafeCFL;
};


template <class VectorType, int dim>
class F77ELCGFMBoundary : public F77GFMBoundary<VectorType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  typedef F77GFMBoundary<VectorType,dim> base;
public:
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::grid_data_type grid_data_type;
  typedef typename base::generic_func_type generic_func_type;
  typedef typename base::point_type point_type;
  typedef AMRELCGFMSolver<VectorType,FixupType,FlagType,dim> amr_elc_solver;

  F77ELCGFMBoundary(generic_func_type bnd, generic_func_type trs, amr_elc_solver& solver) : 
    base(bnd,trs), _solver(solver) {
    base::SetNAux(base::Dim());
  }

  virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi, const BBox& bb,
                           const VectorType* u, DataType* aux, const int& Level, 
                           double t, const int& nc, const int* idx, 
                           const point_type* xc, const DataType* distance, 
                           const point_type* normal) {
    _solver.ExtractBoundaryVelocities(gdu.bbox(),nc,xc,reinterpret_cast<point_type *>(aux));
  }

protected:
  amr_elc_solver& _solver;
};


#endif

---