fsi/sfc-amroc/TubeCJBurnSlot/src/Problem.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROCFLUIDPROBLEMH
#define AMROCFLUIDPROBLEMH

#include "eulerznd3.h"

#undef  NAUX
#define NAUX  2

#include "ClpProblem.h"

#define MaxIntPoints (20)
#define MaxDPoints (10)

#define OWNFLAGGING
#define OWNGFMAMRSOLVER
#include "ClpStdGFMProblem.h" 
#include "AMRGFMInterpolation.h"
#include "StatCPTLevelSet.h"

class FlaggingSpecific : 
  public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
  typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
public:
  FlaggingSpecific(base::solvertype& 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 fflgout
      base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(fflgout));
      base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(fflgout));
      base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(fflgout));
      base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,fflgout));
      base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,fflgout));
#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 registerat(ControlDevice& Ctrl, const std::string& prefix) {
    base::registerat(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 registerat(ControlDevice& Ctrl) {
    registerat(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);
        endforall    
      }
  }

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


template <class DataType, int dim>
class F77CPTLevelSet : public StatCPTLevelSet<DataType,dim> {
  typedef StatCPTLevelSet<DataType,dim> base;
public:
  typedef typename base::gridfcttype gridfcttype;  
  typedef typename base::griddatatype griddatatype;
  typedef genericfortranfunc genericfunctype; 

  typedef void (*lset3functype) ( 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(), flset(0) {}
  F77CPTLevelSet(genericfunctype lset) : base(), flset(lset) {}

  virtual void SetPhi(gridfcttype& 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);    
    endforall
  }

  virtual void SetGrid(griddatatype& gdphi, const int& level, const double& t) {   
    assert(flset != (genericfunctype) 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);
    }

    ((lset3functype) flset)(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(genericfunctype lset) { flset = lset; }
  genericfunctype GetFunc() const { return flset; }

protected:
  genericfunctype flset;
};



class SolverSpecific : 
  public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolationtype;
public:
  typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> outputtype;
  typedef interpolationtype::pointtype pointtype;

  SolverSpecific(IntegratorSpecific& integ, 
                 base::initialconditiontype& init,
                 base::boundaryconditionstype& bc) : base(integ, init, bc) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(fprolong, frestrict));
    SetFileOutput(new outputtype(*this,fflgout)); 
    SetFixup(new FixupSpecific(IntegratorSpecific));
    SetFlagging(new FlaggingSpecific(*this)); 
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(fibndrfl,fitrans),
              new F77CPTLevelSet<DataType,DIM>(flset)));
    Interpolation = new interpolationtype(*this);
    NIntPoints = 0;
    IntName = "ptrack.txt";
  }  
 
  ~SolverSpecific() {
    DeleteGFM(GFM[0]);
    delete Flagging;
    delete Fixup;
    delete Interpolation;
  }

  virtual void registerat(ControlDevice& Ctrl, const std::string& prefix) {
    base::registerat(Ctrl,prefix);
    IntCtrl = base::LocCtrl.getSubDevice("TrackPressure");
    RegisterAt(IntCtrl,"NPoints",NIntPoints);
    char VariableName[32];
    for (int nc=0; nc<MaxIntPoints; nc++) {
      for (int d=0; d<base::Dim(); d++) {
        std::sprintf(VariableName,"Point(%d,%d)",nc+1,d+1);
        RegisterAt(IntCtrl,VariableName,IntPoints[nc](d)); 
        IntPoints[nc](d) = 0.0;
      }
    }
    RegisterAt(IntCtrl,"FileName",IntName);
  } 
  virtual void registerat(ControlDevice& Ctrl) {
    base::registerat(Ctrl);
  }

  virtual void SetupData() {
    base::SetupData();
    Interpolation->SetupData(base::PGH(), base::NGhosts());
  }

  virtual void Advance(double& t, double& dt) {
    base::Advance(t,dt);
    if (NIntPoints<=0 || NIntPoints>MaxIntPoints) return;

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

    DataType* p = new DataType[NIntPoints];
    Interpolation->PointsValues(base::Work(),base::t[0],NIntPoints,IntPoints,p);
    Interpolation->ArrayCombine(VizServer,NIntPoints,p);

    // Append to output file only on processor VizServer
    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      outfile.open(IntName.cstr(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << base::t[0];
      for (int ns=0; ns<NIntPoints; ns++)
        *out << " " << p[ns]; 
      *out << std::endl;
      outfile.close();
      delete out;      
    }

    delete [] p;
  } 

protected:
  IntegratorSpecific IntegratorSpecific;
  InitialConditionSpecific InitialConditionSpecific;
  BoundaryConditionsSpecific BoundaryConditionsSpecific;
  interpolationtype* Interpolation; 
  ControlDevice IntCtrl;
  int NIntPoints;
  pointtype IntPoints[MaxIntPoints];