weno/applications/euler/3d/Homogeneous/src/Problem.h

// -*- C++ -*-

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#define DIM  3
#define NEQUATIONS 9  // Euler equations for gases in 3D (5 fields), 
                       // +1 fields for passive scalars
                       // +1 field for temperature
                       // +1 feild to output dcflag and 
                       // +1 fieldsgske
#define NVARS       6
  
#include "WENOProblem.h"
#define OWN_AMRSOLVER
#include "WENOStdProblem.h"

class SolverSpecific : 
  public AMRSolver<VectorType,FixupType,FlagType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef WENOF77FileOutput<VectorType,DIM> output_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) :
    AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
#ifdef f_flgout
    SetFileOutput(new WENOF77FileOutput<VectorType,DIM>(f_flgout,f_bounds)); 
#else   
    SetFileOutput(new FileOutput<VectorType,DIM>()); 
#endif
    SetFixup(new FixupSpecific());
    SetFlagging(new FlaggingSpecific(*this)); 
  }  

  ~SolverSpecific() {
    delete _LevelTransfer;
    delete _Flagging;
    delete _Fixup;
    delete _FileOutput;
  }

  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);
  }

  virtual void SetupData() {
    base::SetupData();
  }

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

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

    // statistics 
    int irho = 0;
    int iu   = 1;
    int iv   = 2;
    int iw   = 3;
    int isgs = 11;
    int ip = 7;
    int Level = 0 ;
    int Time = CurrentTime(base::GH(),Level);
    int TStep = TimeStep(base::U(),Level);
    double rkin, skin, pavg, pvar, rhoavg, rhovar, vol;

    vol = base::GH().wholebbox().upper(0)-base::GH().wholebbox().lower(0)+1;
    vol *= base::GH().wholebbox().upper(1)-base::GH().wholebbox().lower(1)+1;
    vol *= base::GH().wholebbox().upper(2)-base::GH().wholebbox().lower(2)+1;

    // subgrid kinetic energy
    forall (base::U(),Time,Level,c)
      ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), 
                                                    Time, Level, isgs, 
                                                    base::t[Level]);
    end_forall
    skin = Work().GF_sum(Time,Level)/vol;

    // pressure
    forall (base::U(),Time,Level,c)
      ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), 
                                                    Time, Level, ip, 
                                                    base::t[Level]);
       Coords ss(base::U()(Time,Level,c).bbox().stepsize());
       BBox bbi(base::U().interiorbbox(Time,Level,c));
       BeginFastIndex3(Work, Work()(Time,Level,c).bbox(), 
                       Work()(Time,Level,c).data(), DataType);
       BeginFastIndex3(WorkN, Work()(Time+TStep,Level,c).bbox(), 
                        Work()(Time+TStep,Level,c).data(), DataType);
       DataType p;
       for_3 (i, j, k, bbi, ss)
         p = FastIndex3(Work,i,j,k);
       FastIndex3(WorkN,i,j,k) = p*p;
       end_for
       EndFastIndex3(Work);
       EndFastIndex3(WorkN);
    end_forall
    pavg = Work().GF_sum(Time,Level)/vol;
    pvar = Work().GF_sum(Time+TStep,Level)/vol-pavg*pavg;

    // density
    forall (base::U(),Time,Level,c)
       Coords ss(base::U()(Time,Level,c).bbox().stepsize());
       BBox bbi(base::U().interiorbbox(Time,Level,c));
       BeginFastIndex3(U, base::U()(Time,Level,c).bbox(), 
                       base::U()(Time,Level,c).data(), VectorType);
       BeginFastIndex3(Work, Work()(Time,Level,c).bbox(), 
                       Work()(Time,Level,c).data(), DataType);
       BeginFastIndex3(WorkN, Work()(Time+TStep,Level,c).bbox(), 
                        Work()(Time+TStep,Level,c).data(), DataType);
       DataType rho;
       for_3 (i, j, k, bbi, ss)
         rho = FastIndex3(U,i,j,k)(irho);
         FastIndex3(Work,i,j,k) = rho;
         FastIndex3(WorkN,i,j,k) = rho*rho;
       end_for
       EndFastIndex3(U);
       EndFastIndex3(Work);
       EndFastIndex3(WorkN);
    end_forall
    rhoavg = Work().GF_sum(Time,Level)/vol;
    rhovar = Work().GF_sum(Time+TStep,Level)/vol-rhoavg*rhoavg;

    // reslved kinetic energy
    forall (base::U(),Time,Level,c)
       Coords ss(base::U()(Time,Level,c).bbox().stepsize());
       BBox bbi(base::U().interiorbbox(Time,Level,c));
       BeginFastIndex3(U, base::U()(Time,Level,c).bbox(), 
                       base::U()(Time,Level,c).data(), VectorType);
       BeginFastIndex3(Work, Work()(Time,Level,c).bbox(), 
                       Work()(Time,Level,c).data(), DataType);
       DataType rho, u, v, w;
       for_3 (i, j, k, bbi, ss)
         rho = FastIndex3(U,i,j,k)(irho);
         u = FastIndex3(U,i,j,k)(iu)/rho;
         v = FastIndex3(U,i,j,k)(iv)/rho;
         w = FastIndex3(U,i,j,k)(iw)/rho;
         FastIndex3(Work,i,j,k) = (u*u+v*v+w*w)/2.0;
       end_for
       EndFastIndex3(U);      
       EndFastIndex3(Work);
    end_forall
    rkin = Work().GF_sum(Time,Level)/vol;

    if (MY_PROC == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "stats.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << base::t[0] << " " << rkin << " " << skin << " " 
           << (rkin+skin) << " " << pvar << " " << rhovar << std::endl;
      outfile.close();
      delete out;      
    }
         
    return cfl;
  } 

protected: