clawpack/applications/euler/3d/ConvShock_Cyl/src/Problem.h

// -*- C++ -*-

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

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include "eulerznd3.h"
#include "ClpProblem.h"

#define f_track FORTRAN_NAME(track, TRACK)
#define f_lsetrfl FORTRAN_NAME(lsrfl, LSRFL)
#define f_lsetrdi FORTRAN_NAME(lsrdi, LSRDI)

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

#define OWN_GFMAMRSOLVER
#include "ClpStdGFMProblem.h"
#include "AMRGFMInterpolation.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 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_ibndex,f_itrans),
              new F77GFMLevelSet<DataType,DIM>(f_lset)));
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
              new F77GFMLevelSet<DataType,DIM>(f_lsetrdi)));
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
              new F77GFMLevelSet<DataType,DIM>(f_lsetrfl)));
    _Interpolation = new interpolation_type(*this);
  }  
 
  ~SolverSpecific() {
    DeleteGFM(_GFM[0]);
    DeleteGFM(_GFM[1]);
    DeleteGFM(_GFM[2]);
    delete _Flagging;
    delete _Fixup;
    delete _Interpolation;
  }

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

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

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

//     int me = MY_PROC; 
//     int Npoints = 100, npos;
//     double pi = 4.0*std::atan(1.0), alpha = 0.0, pmax = 0.0;
//     double shk_pos[10];
//     double fmin = 0.1, fmax = 1.e4, dfmin = 10.;
//     for (npos=0; npos<10; alpha+=5.0, npos++) {

//       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] = 8.0/Npoints*n;
//      xc[n](0) = r[n]*std::cos(alpha*pi/180.); xc[n](1) = r[n]*std::sin(alpha*pi/180.);
//      xc[n](2) = 0.0;
//       }
    
//       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]);
//       }

//       _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];
//      if (base::t[0]<0.28) {
//        dfmin = 10.;
//        f_track(Npoints,p,r,xr,nc,fmin,fmax,dfmin);
//        shk_pos[npos] = xr[0];
//      }
//      else {
//        dfmin = 1000.;
//        f_track(Npoints,p,r,xr,nc,fmin,fmax,dfmin);
//        shk_pos[npos] = xr[nc-1];
//      }
//       }
      
//       delete [] p;
//       delete [] r;
//       delete [] xr;  
//       delete [] xc;
//     }

//     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());
//       double sum = 0.0;
//       for (npos=0; npos<10; npos++) 
//      sum += shk_pos[npos];
//       *out << base::t[0] << " " << pmax << " " << sum/10.;      
//       for (npos=0; npos<10; npos++) 
//      *out << " " << shk_pos[npos];
//       *out << std::endl;
//       outfile.close();
//       delete out;      
//     }

//     return cfl;
  } 

protected:
  interpolation_type* _Interpolation;