fsi/sfc-amroc/WaterBlastFracture/src/FluidProblem.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_FLUIDPROBLEM_H
#define AMROC_FLUIDPROBLEM_H

#include "eulerm3.h"

#include "ClpProblem.h"

#define f_combl_read FORTRAN_NAME(cblread, CBLREAD)
#define f_combl_write FORTRAN_NAME(cblwrite, CBLWRITE)
#define f_ipvel_piston FORTRAN_NAME(ipvelpiston, IPVELPISTON)
#define f_lset_piston FORTRAN_NAME(lspiston, LSPISTON)

extern "C" {
  void f_combl_read(DOUBLE& xm, DOUBLE& vm, DOUBLE& cm, DOUBLE& rd, INTEGER& Np); 
  void f_combl_write(DOUBLE& xm, DOUBLE& vm, DOUBLE& cm, DOUBLE& rd, INTEGER& Np); 
  void f_ipvel_piston(); 
  void f_lset_piston(); 
}

#define MaxDPoints (10)
#define OUTPUTPOINTSMAX  20

#define OWN_FLAGGING
#define OWN_ELCGFMAMRSOLVER
#include "ClpStdELCGFMProblem.h" 
#include "AMRGFMInterpolation.h"


class FlaggingSpecific : 
  public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
  typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
public:
  FlaggingSpecific(base::solver_type& solver) : base(solver) {
      base::AddCriterion(new ByValue<VectorType,FlagType,DIM>());
      base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
      base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
      base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
      base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
#ifdef f_flgout
      base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
      base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
#endif
      dcmin = new int[MaxDPoints*base::Dim()];
      dcmax = new int[MaxDPoints*base::Dim()];
      dclev = new int[MaxDPoints];
      for (int d=0; d<MaxDPoints*base::Dim(); d++) {
        dcmin[d] = 0; dcmax[d] = -1;
      }
      for (int i=0; i<MaxDPoints; i++)
        dclev[i] = 1000000;
    }

  ~FlaggingSpecific() { 
    DeleteAllCriterions(); 
    delete [] dcmin;
    delete [] dcmax;
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl, prefix);
    char VariableName[32];
    for (int i=0; i<MaxDPoints; i++) {
      for (int d=0; d<base::Dim(); d++) {
        sprintf(VariableName,"DCellsMin(%d,%d)",i+1,d+1);
        RegisterAt(base::LocCtrl,VariableName,dcmin[i*base::Dim()+d]);
        sprintf(VariableName,"DCellsMax(%d,%d)",i+1,d+1);
        RegisterAt(base::LocCtrl,VariableName,dcmax[i*base::Dim()+d]);
      }
      sprintf(VariableName,"DMinLevel(%d)",i+1);
      RegisterAt(base::LocCtrl,VariableName,dclev[i]);
    }
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual void SetFlags(const int Time, const int Level, double t, double dt) {
    base::SetFlags(Time, Level, t, dt);
    for (int i=0; i<MaxDPoints; i++) 
      if (Level>=dclev[i]) {
        Coords lb(base::Dim(), dcmin+i*base::Dim()), ub(base::Dim(), dcmax+i*base::Dim()), 
                  s(base::Dim(),1);
        BBox bbox(lb*RefinedBy(base::GH(),MaxLevel(base::GH())),
                  ub*RefinedBy(base::GH(),MaxLevel(base::GH())),
                  s*StepSize(base::GH(),Level));
        forall (base::Flags(),Time,Level,c)  
          base::Flags()(Time,Level,c).equals(GoodPoint, bbox);
        end_forall    
      }
  }

private:
  int *dcmin, *dcmax, *dclev;
};


template <class DataType, int dim>
class F77CPTLevelSet : public CPTLevelSet<DataType,dim> {
  typedef CPTLevelSet<DataType,dim> base;
public:
  typedef typename base::grid_fct_type grid_fct_type;  
  typedef typename base::grid_data_type grid_data_type;
  typedef generic_fortran_func generic_func_type; 

  typedef void (*lset_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
                                     const INTEGER& mbc,
                                     const INTEGER& mx, const INTEGER& my, const INTEGER& mz,
                                     const DOUBLE x[], const DOUBLE y[], const DOUBLE z[],
                                     const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz, 
                                     DataType q[], const DOUBLE& t);

  F77CPTLevelSet() : base(), f_lset(0) {}
  F77CPTLevelSet(generic_func_type lset) : base(), f_lset(lset) {}

  virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) { 
    base::SetPhi(phi,Time,Level,t);
    forall (phi,Time,Level,c)
      SetGrid(phi(Time,Level,c),Level,t);    
    end_forall
  }

  virtual void SetGrid(grid_data_type& gdphi, const int& level, const double& t) {   
    assert(f_lset != (generic_func_type) 0);
    Coords ex = gdphi.extents();
    DCoords lbcorner = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
    DCoords dx = base::GH().worldStep(gdphi.stepsize());
    int maxm[3], mx[3], d;
    DataType* x[3];
    for (d=0; d<dim; d++) {
      maxm[d] = ex(d);
      mx[d] = ex(d)-2*base::NGhosts();
      x[d] = new DataType[maxm[d]];
      for (int i=0; i<maxm[d]; i++) x[d][i] = (i+0.5)*dx(d)+lbcorner(d);
    }

    ((lset_3_func_type) f_lset)(AA(3,mx), base::NGhosts(), AA(3,mx), 
                                AA(3,x), AA(3,dx), gdphi.data(), t);

    for (d=0; d<dim; d++) 
      delete [] x[d];
  }

  inline void SetFunc(generic_func_type lset) { f_lset = lset; }
  generic_func_type GetFunc() const { return f_lset; }

protected:
  generic_func_type f_lset;
};


class FluidSolverSpecific : 
  public AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
public:
  typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
  typedef base::point_type point_type;
  typedef base::eul_comm_type eul_comm_type;

  FluidSolverSpecific() : base(_IntegratorSpecific, _InitialConditionSpecific, 
                               _BoundaryConditionsSpecific) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
    SetFileOutput(new output_type(*this,f_flgout)); 
    SetFixup(new FixupSpecific(_IntegratorSpecific));
    SetFlagging(new FlaggingSpecific(*this)); 
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77ELCGFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,*this),
              new CPTLevelSet<DataType,DIM>()));
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
              new F77GFMLevelSet<DataType,DIM>(f_lset)));
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,f_ipvel_piston),
              new F77GFMLevelSet<DataType,DIM>(f_lset_piston)));
    SetCoupleGFM(0);
    _Interpolation = new interpolation_type(*this);
  }  
 
  ~FluidSolverSpecific() {
    DeleteGFM(_GFM[2]);
    DeleteGFM(_GFM[1]);
    DeleteGFM(_GFM[0]);
    delete _Flagging;
    delete _Fixup;
    delete _Interpolation;
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
    SensorsCtrl = base::LocCtrl.getSubDevice("SensorsOutput");
    RegisterAt(SensorsCtrl,"OutputEvery",OutputEvery);
    RegisterAt(SensorsCtrl,"NSensors",Nxc);
    char VariableName[32];
    for (register int i=0; i<OUTPUTPOINTSMAX; i++) 
      for (register int d=0; d<base::Dim(); d++) {
        std::sprintf(VariableName,"xc(%d,%d)",i+1,d+1);
        RegisterAt(SensorsCtrl,VariableName,xc[i](d));
      }
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    base::register_at(Ctrl);
  }

  virtual void SetupData() {
    base::SetupData();
    _Interpolation->SetupData(base::PGH(), base::NGhosts());
    if (OutputEvery<0) OutputEvery=0;
    if (Nxc<0) Nxc=0;
    if (Nxc>OUTPUTPOINTSMAX) Nxc=OUTPUTPOINTSMAX;
  }

  virtual void Restart(double& t, double& dt) {
    std::ifstream ifs("boundary.cp", std::ios::in);
    int Np;
    double xm, vm, cm, rd;
    f_combl_read(xm,vm,cm,rd,Np);
    ifs >> xm >> vm;
    ifs.close();
    f_combl_write(xm,vm,cm,rd,Np);
    base::Restart(t,dt);
  }
  virtual void Checkpointing() {
    int me = MY_PROC; 
    if (me == VizServer) {
      int Np;
      double xm, vm, cm, rd;
      f_combl_read(xm,vm,cm,rd,Np);
      std::ofstream ofs("boundary.cp", std::ios::out);
      ofs << xm << " " << vm << std::endl; 
      ofs.close();
    }
    base::Checkpointing();
  }

  virtual void SendBoundaryData() {
    START_WATCH
      for (register int l=0; l<=FineLevel(base::GH()); l++) 
        ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), 
                                                CurrentTime(base::GH(),l), l, 
                                                base::Dim()+4, base::t[l]);
    END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
    base::SendBoundaryData();
  }

  virtual void ModifyELCSendData(eul_comm_type* eulComm) {
    DataType distance = -GFM(CoupleGFM()).Boundary().Cutoff();
    int Npoints = _eulComm->getNumberOfFaces();
    double *press_data = _eulComm->getPressuresData();
    const int* con = _eulComm->getConnectivitiesData();
    const point_type *pos = reinterpret_cast<const point_type *>(_eulComm->getPositionsData());
    for (register int n=0; n<Npoints; n++) {
      bool inside = true;
      for (register int d=0; d<base::Dim(); d++) {
        const point_type &p = pos[con[base::Dim()*n+d]]; 
        if (std::sqrt(p(1)*p(1)+p(2)*p(2))>0.0323) {
          inside = false;
          break;
        }
      }
      if (!inside) 
        press_data[n] = std::numeric_limits<DataType>::max();
    }
  }
  
  // Overloading of Advance() to add pointwise output after each level 0 time step
  // and propagate piston
  virtual void Advance(double& t, double& dt) {

    int me = MY_PROC; 
    DataType force_old;
    if (base::GFM(2).IsUsed())
      CalculateForce(force_old);

    base::Advance(t,dt);

    if (base::GFM(2).IsUsed()) {
      DataType a, force_new;
      CalculateForce(force_new);
      
      int Np;
      double xm, vm, cm, rd;
      f_combl_read(xm,vm,cm,rd,Np);
      
      a = 0.5/cm*(force_old+force_new);
      vm += a*dt;
      xm += vm*dt;
      
      f_combl_write(xm,vm,cm,rd,Np);
    
      // Append to output file only on processor VizServer
      if (me == VizServer) {
        std::ofstream outfile;
        std::ostream* out;
        std::string fname = "boundary.txt";
        outfile.open(fname.c_str(), std::ios::out | std::ios::app);
        out = new std::ostream(outfile.rdbuf());
        *out << t << " " << force_new << " " << xm << " " << vm << std::endl;
        outfile.close();
        delete out;      
      }
    }

    if (!OutputEvery) return;

    int TimeZero  = CurrentTime(base::GH(),0);
    TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
    if (TimeZero % OutputEvery != 0) return;

    VectorType* dat = new VectorType[Nxc];

    // Interpolate/extrapolate pressure values and assemble results on processor VizServer
    _Interpolation->PointsValues(base::U(),t,Nxc,xc,dat);
    int NValues = Nxc*base::NEquations();
    _Interpolation->ArrayCombine(VizServer,NValues,&(dat[0](0)));

    // Append to output file only on processor VizServer
    if (MY_PROC == VizServer) {
      DataType* output_data = new DataType[_FileOutput->Ncnt()*Nxc];
      BBox bb(3,0,0,0,Nxc-1,0,0,1,1,1);
      vec_grid_data_type val(bb,dat);
      for (int cnt=1;cnt<=_FileOutput->Ncnt(); cnt++) {
        if (_FileOutput->OutputName(cnt-1).c_str()[0] == '-') continue;
        grid_data_type output(bb,output_data+(cnt-1)*Nxc);
        ((output_type::out_3_func_type) (((output_type*) base::_FileOutput)->GetOutFunc()))
          (FA(3,val),FA(3,output),BOUNDING_BOX(bb),base::NEquations(),cnt,t);
        DataType* dummy;
        output.deallocate(dummy);
      }
      VectorType* dummy;
      val.deallocate(dummy);
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "sensors.txt";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << t << " ";
      for (register int j=0;j<Nxc; j++)
        for (int cnt=0;cnt<_FileOutput->Ncnt(); cnt++) 
          if (_FileOutput->OutputName(cnt).c_str()[0] != '-') 
            *out << " " << output_data[cnt*Nxc+j];
      *out << std::endl;
      delete [] output_data;
      outfile.close();
      delete out;      
    }
    delete [] dat;
  } 

  void CalculateForce(DataType& sum) {
    int me = MY_PROC; 

    double pi = 4.0*std::atan(1.0);
    
    // Read sphere info from F77 common block
    int Np;
    double xm, vm, cm, rd;
    f_combl_read(xm,vm,cm,rd,Np);

    // Calculate pressure
    for (register int l=0; l<=FineLevel(base::GH()); l++) {
      int Time = CurrentTime(base::GH(),l);
      int press_cnt = base::Dim()+4;      
      ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
                                              press_cnt, base::t[l]);
    }

    int Npoints = Np*Np;
    DataType* p = new DataType[Npoints];
    point_type* xc = new point_type[Npoints];

    register int n;

    for (n=0; n<Np; n++) 
      for (register int m=0; m<Np; m++) {
        xc[n*Np+m](0) = xm; 
        xc[n*Np+m](1) = (2.*n/Np-1.)*rd;
        xc[n*Np+m](2) = (2.*m/Np-1.)*rd;
      }
    
    // Interpolate/extrapolate pressure values and assemble results on processor VizServer
    _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
    _Interpolation->ArrayCombine(VizServer,Npoints,p);

    // Do the integration on the surface only on processor VizServer
    if (me == VizServer) {
      int vp = 0;
      sum = 0.;
      for (n=0; n<Npoints; n++)
        if (p[n]>-1.e37) {
          sum += p[n]; vp++;
        }
      sum /= vp;
      sum -= 101325.;
      sum *= pi*rd*rd;

      delete [] p;
      delete [] xc;
    }

#ifdef DAGH_NO_MPI
#else
    if (comm_service::dce() && comm_service::proc_num() > 1) {  
      int num = comm_service::proc_num();
      MPI_Status status;
      int R;
      int size = sizeof(DataType);
      int tag = 10001;
      if (me == VizServer) 
        for (int proc=0; proc<num; proc++) {
          if (proc==VizServer) continue;
          R = MPI_Send((void *)&sum, size, MPI_BYTE, proc, tag, comm_service::comm());
          if ( MPI_SUCCESS != R ) 
            comm_service::error_die( "SumCombine", "MPI_Send", R );
        }
      else {
        R = MPI_Recv((void *)&sum, size, MPI_BYTE, VizServer, tag, comm_service::comm(), &status);
        if ( MPI_SUCCESS != R ) 
          comm_service::error_die( "ArrayCombine", "MPI_Recv", R );
      }
    }
#endif
  }

protected:
  IntegratorSpecific _IntegratorSpecific;
  InitialConditionSpecific _InitialConditionSpecific;
  BoundaryConditionsSpecific _BoundaryConditionsSpecific;
  int OutputEvery, Nxc;
  point_type xc[OUTPUTPOINTSMAX];
  interpolation_type* _Interpolation;