weno/applications/euler/3d/RM_AirSF6/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 feild for temperature
                       // +1 feild to output dcflag and 
                       // +1 fieldsgske
#define NVARS     6   // Fixup used only for vector of state and 2 scalars
#define OWN_STATISTICS
#define OWN_AMRSOLVER
#include "WENOProblem.h"
#include "WENOStdProblem.h"
#ifdef OWN_STATISTICS
#include "AMRInterpolation.h"
#include "WENOStatistics.h"
#endif

#ifdef OWN_AMRSOLVER

#define f_init_commongm FORTRAN_NAME(comblgm, COMBLGM)

extern "C" {
  void f_init_commongm(const INTEGER& meqn, INTEGER* shape, DOUBLE* geom, 
                       INTEGER& maxlev, INTEGER* refine, INTEGER* per);
}

class SolverSpecific : 
  public AMRSolver<VectorType,FixupType,FlagType,DIM> {
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef WENOFixup<VectorType,FixupType,DIM> weno_fixup_type;
  typedef WENOIntegrator<VectorType,DIM> weno_integ_type;
#ifdef OWN_STATISTICS
  typedef AMRInterpolation<VectorType,DIM> interpolation_type;
  typedef WENOF77FileOutput<VectorType,DIM> output_type;
  typedef WENOStatistics<VectorType,interpolation_type,output_type,DIM> stat_type;
#endif
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)); 
#ifdef OWN_STATISTICS
    _Interpolation = new interpolation_type();
    _Stats = new stat_type(_Interpolation, (output_type*)_FileOutput);
#endif
  }

#ifdef OWN_STATISTICS
    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);
        if (_Stats) _Stats->register_at(base::LocCtrl, ""); 
    }
#endif
    
  virtual void SetupData() {
    base::SetupData();
    int dim = DIM;
    int d, per[DIM];
    for (d=0; d<dim; d++) { 
      per[d] = base::GH().periodicboundary(d);
    }
    f_init_commongm(dim,shape,geom,MaxLev,_RefineFactor,per);
#ifdef OWN_STATISTICS
    _Interpolation->SetupData(base::PGH(), base::NGhosts());
    _Stats->Setup(base::PGH(), base::NGhosts(), base::shape, base::geom);
#endif
  }

#ifdef OWN_STATISTICS
    virtual void Advance(double& t, double& dt) {
        base::Advance(t, dt);
        _Stats->Evaluate(base::t, base::U(), base::Work());
    }

   virtual void Initialize_(const double& dt_start) { 
      base::Initialize_(dt_start);
      _Stats->Evaluate(base::t, base::U(), base::Work());
  }
#endif

  ~SolverSpecific() {
    delete _LevelTransfer;
    delete _Flagging;
    delete _Fixup;
    delete _FileOutput;
#ifdef OWN_STATISTICS
    delete _Interpolation;
    delete _Stats;
#endif
  }
#ifdef OWN_STATISTICS
protected:
  interpolation_type* _Interpolation; 
  stat_type* _Stats;
#endif
};