weno/applications/euler/3d/Jet/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   // Fixup used only for vector of state and 2 scalars

//#define SYNCTIME
#define OWN_AMRSOLVER
#define OWN_FLAGGING

#include "WENOProblem.h"
#ifdef SYNCTIME
#include "WENOStdSyncTimeProblem.h"
#else
#include "WENOStdProblem.h"
#endif
#include "AMRInterpolation.h"
#include "WENOStatistics.h"

#ifdef OWN_FLAGGING
#define NZONES 2

class FlaggingSpecific : 
  public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
  typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
public:
  FlaggingSpecific(base::solver_type& 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 f_flgout
      base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
      base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
#endif
    }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl, prefix);
    char VariableName[32];
    for (int iz=0; iz < NZONES ; iz++)
      for (int d=0; d<DIM; d++) {
        sprintf(VariableName,"DCellsMin%d(%d)",iz+1,d+1);
        RegisterAt(base::LocCtrl,VariableName,dcmin[iz][d]);
        sprintf(VariableName,"DCellsMax%d(%d)",iz+1,d+1);
        RegisterAt(base::LocCtrl,VariableName,dcmax[iz][d]);
      } 
    for (int d=0; d<DIM; d++) {
      sprintf(VariableName,"FCellsMin(%d)",d+1);
      RegisterAt(base::LocCtrl,VariableName,scmin[d]);
      sprintf(VariableName,"FCellsMax(%d)",d+1);
      RegisterAt(base::LocCtrl,VariableName,scmax[d]);
    } 
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual void SetFlags(const int Time, const int Level, double t, double dt) {
    base::SetFlags(Time, Level, t, dt);
    int max_level = MaxLevel(base::GH());
    for ( int iz = 0 ; iz < NZONES ; iz ++ ) {
      Coords lb(base::Dim(), dcmin[iz]), ub(base::Dim(), dcmax[iz]), 
        s(base::Dim(),1);
      BBox bbox(lb*StepSize(base::GH(),0), ub*StepSize(base::GH(),0),
                s*StepSize(base::GH(),0));
      bbox.refine(StepSize(base::GH(),0)/StepSize(base::GH(),Level));
      
      forall (base::Flags(),Time,Level,c)  
        if ( Level > max_level - (2+iz) )
          base::Flags()(Time,Level,c).equals(GoodPoint, bbox);
      end_forall
    }
    // Mark the slot
    Coords lb(base::Dim(), scmin), ub(base::Dim(), scmax), s(base::Dim(),1);
    BBox bbox(lb*StepSize(base::GH(),0), ub*StepSize(base::GH(),0),
              s*StepSize(base::GH(),0));
    bbox.refine(StepSize(base::GH(),0)/StepSize(base::GH(),Level));
      
    forall (base::Flags(),Time,Level,c)
      base::Flags()(Time,Level,c).equals(BadPoint, bbox);
    end_forall
  }

  ~FlaggingSpecific() { DeleteAllCriterions(); }
private:
  int dcmin[NZONES][DIM], dcmax[NZONES][DIM];
  int scmin[DIM], scmax[DIM];
};
#endif

#ifdef OWN_AMRSOLVER

class SolverSpecific : 
#ifdef SYNCTIME
  public AMRSolverSyncTime<VectorType,FixupType,FlagType,DIM> {
  typedef AMRSolverSyncTime<VectorType,FixupType,FlagType,DIM> base;
#else
  public AMRSolver<VectorType,FixupType,FlagType,DIM> {
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
#endif
  typedef AMRInterpolation<VectorType,DIM> interpolation_type;
  typedef WENOF77FileOutput<VectorType,DIM> output_type;
  typedef WENOStatistics<VectorType,interpolation_type,output_type,DIM> stat_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) :
#ifdef SYNCTIME
    AMRSolverSyncTime<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
#else
    AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
#endif
    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)); 
    _Interpolation = new interpolation_type();
    _Stats = new stat_type(_Interpolation, (output_type*)_FileOutput);
  }  

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

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

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

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

protected:
  interpolation_type* _Interpolation; 
  stat_type* _Stats;