clawpack/applications/euler/2d/Conical_Shocktube/src/Problem.h

// -*- C++ -*-

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include "euler2.h"

#undef  NAUX
#define NAUX       1

#include "ClpProblem.h"

#define f_track FORTRAN_NAME(track, TRACK)

extern "C" {
  void f_track(const INTEGER &, const DOUBLE *, const DOUBLE *, DOUBLE *, INTEGER&,
               const DOUBLE&, const DOUBLE&, const DOUBLE&,
               DOUBLE&, DOUBLE&, DOUBLE&, DOUBLE&, DOUBLE&, DOUBLE&); 
}

#define OWN_GFMAMRSOLVER
#include "ClpStdGFMProblem.h"
#include "AMRGFMInterpolation.h"
#include "SpecialOutput/AMRGridAccumulation.h"

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> interpolation_type;
  typedef AMRGridAccumulation<DataType,DIM> accumulation_type;
  typedef interpolation_type::point_type point_type;
  typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) : base(integ, init, bc) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
#ifdef f_flgout
    SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout)); 
#else   
    SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this)); 
#endif
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
              new F77GFMLevelSet<DataType,DIM>(f_lset)));
    _Interpolation = new interpolation_type(*this);
    _Accumulation = new accumulation_type();
    MaxType = 1;
  }  
 
  ~SolverSpecific() {
    DeleteGFM(_GFM[0]);
    delete _Flagging;
    delete _Fixup;
    delete _Interpolation;
    delete _Accumulation;
  }

  virtual void register_at(ControlDevice& Ctrl) { base::register_at(Ctrl); }
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl, prefix);
    RegisterAt(base::LocCtrl,"MaxType",MaxType);
    _Accumulation->register_at(base::LocCtrl, prefix);
  }

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

  virtual void Initialize_(const double& dt_start) {
    base::Initialize_(dt_start);
    _Accumulation->Initialize();
  }

  virtual void Output() {
    int Time = CurrentTime(base::GH(),0);
    if (Time == base::LastOutputTime()) return;
    _Accumulation->GridCombine(VizServer);
    int me = MY_PROC; 
    if (me == VizServer) {
      int Level = _Accumulation->AccLevel();
      base::FileOutput_().WriteOut(_Accumulation->AccGrid(),_Accumulation->AccName(),
                                   Time,Level,base::t[0]);
    }
    base::Output();
  }

  virtual void Checkpointing_(const char* CheckpointFile) {
    base::Checkpointing_(CheckpointFile);
    _Accumulation->Checkpointing();
  }

  virtual bool Restart_(const char* CheckpointFile) {
    bool ret = base::Restart_(CheckpointFile);
    if (ret) _Accumulation->Restart();
    return ret;
  }

  virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone, 
                            bool ShadowAllowed, bool DoFixup,
                            bool RecomposeBaseLev, bool RecomposeHighLev) {

    base::AdvanceLevel(Level,RegridEvery,RegridDone,ShadowAllowed,
                       DoFixup,RecomposeBaseLev,RecomposeHighLev);

    int Time = CurrentTime(base::GH(),Level);
    if (MaxType) {
      AdaptiveBoundaryUpdate(base::U(),Time,Level);
      Sync(base::U(),Time,Level);
      ExternalBoundaryUpdate(base::U(),Time,Level);
      
      int irho = 1;
      int iu   = NEQUSED-DIM-2;
      int iv   = NEQUSED-DIM-1;
      int TStep = TimeStep(base::U(),Level);
      ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time+TStep, 
                                                    Level, irho, base::t[Level]);
      forall (base::U(),Time,Level,c)
        Coords ss(base::U()(Time,Level,c).bbox().stepsize());
        BBox bbi(base::U().interiorbbox(Time,Level,c));
        BeginFastIndex2(U, base::U()(Time,Level,c).bbox(), 
                        base::U()(Time,Level,c).data(), VectorType);
        BeginFastIndex2(Work, Work()(Time,Level,c).bbox(), 
                        Work()(Time,Level,c).data(), DataType);
        BeginFastIndex2(WorkN, Work()(Time+TStep,Level,c).bbox(), 
                        Work()(Time+TStep,Level,c).data(), DataType);
        DataType uy, vx;
        DCoords dx = base::GH().worldStep(ss);
        for_2 (i, j, bbi, ss)
          vx = (FastIndex2(U,i+ss(0),j)(iv)/FastIndex2(WorkN,i+ss(0),j)-
                FastIndex2(U,i-ss(0),j)(iv)/FastIndex2(WorkN,i-ss(0),j))/(2.0*dx(0));
          uy = (FastIndex2(U,i,j+ss(1))(iu)/FastIndex2(WorkN,i,j+ss(1))-
                FastIndex2(U,i,j-ss(1))(iu)/FastIndex2(WorkN,i,j-ss(1)))/(2.0*dx(1));
          FastIndex2(Work,i,j) = std::fabs(vx-uy);
        end_for
        EndFastIndex2(U);      
        EndFastIndex2(Work);      
        EndFastIndex2(WorkN);      
      end_forall
    }
    else {
      int ipress = base::Dim()+4;      
      ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level, 
                                                    ipress, base::t[Level]);
    }
  
    _Accumulation->Accumulate(base::Work(), Time,  Level, base::t[Level]);
  }

  virtual double Tick(int VariableTimeStepping, const double dtv[], const double cflv[],
                      int& Rejections) {

    double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);

    int me = MY_PROC;
    // int Npoints=2000;
    int Npoints; 
    int finelevel = FineLevel(base::GH());
    Npoints = base::shape[0]*RefinedBy(base::GH(),finelevel);
    double dx = DeltaX(base::GH(),0,finelevel);
    //double dy = DeltaX(base::GH(),1,finelevel);

    double pmax = 0.0;
    // double rpmaxloc = -1.0;
    double shk_pos = -0.12;
    double aMa,aMb,aMc,aMd,aMe,aMm;
    double fmin = 150.00, fmax = 20000.;
    double dfmin = 0.1;
    double x0rdi = 0.424779, y0rdi;
    // the probe diameter is .0016meters
    y0rdi = 0.0016;
    // y0rdi = .5*dy;
    double reflection_time = 0.00021;
    // double lngth = 0.424779+0.12;
    // Npoints = int(lngth/dx);
 
    DataType* p = new DataType[Npoints];
    DataType* r = new DataType[Npoints];
    DataType* xr = new DataType[Npoints];
    point_type* xc = new point_type[Npoints];
    register int n;

    for (n=0; n<Npoints; n++) {
      r[n] =  base::geom[0]+Npoints*dx -(n+0.5)*dx;
      //r[n] = lngth/Npoints*n;
      if((r[n])>x0rdi) {
        xc[n](0) = x0rdi;
      } else {
        xc[n](0) = r[n];
      }
      xc[n](1) = y0rdi;
    }

    // get the pressure from the conserved vector of state
    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]);
    }
    // interpolate the pressure values
    _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
    _Interpolation->ArrayCombine(VizServer,Npoints,p);
    
    if (me == VizServer) {
      int nc;
      for (n=0; n<Npoints; n++) 
        if (p[n]>pmax) {
          pmax = p[n];
          // rpmaxloc = x0rdi-r[n];
          // rpmaxloc = xc[n](0);
        }
      if (base::t[0]<reflection_time) {
        dfmin = 0.0000000001;
        f_track(Npoints,p,r,xr,nc,fmin,fmax,dfmin,aMa,aMb,aMc,aMd,aMe,aMm);
        // 0 is the point with the largest x-coordinate (see indexing above)
        // shk_pos = x0rdi-xr[0];
        shk_pos = xr[0];
      }
      else {
        dfmin = 0.1;
        f_track(Npoints,p,r,xr,nc,fmin,fmax,dfmin,aMa,aMb,aMc,aMd,aMe,aMm);
        // nc-1 is the point with the smallest x-coordinate (see indexing above)
        //shk_pos = x0rdi-xr[nc-1];
        shk_pos = xr[nc-1];
      }
    }
      
    delete [] p;
    delete [] r;
    delete [] xr;       
    delete [] xc;
  
    
    // done extracting the shock locations along the ray
    // now output the locations to the file shk.dat

    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "shk.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      // output the time and the max pressure
      *out << base::t[0] << " " << shk_pos << " " << pmax;
      // also output the shock location
      *out << " " << aMa << " " << aMb << " " << aMc;
      *out << " " << aMd << " " << aMe << " " << aMm;
      *out << std::endl;
      outfile.close();
      delete out;      
    }
    return cfl;
  } 

protected:
  int MaxType;
  interpolation_type* _Interpolation; 
  accumulation_type* _Accumulation;