weno/applications/euler_chem/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 the temperature
                       // +1 field to output dcflag and 
                       // +1 field for sgske
#define NVARS     6   // Fixup used only for vector of state and 2 scalars
  
//#define SYNCTIME
//#define OWN_INITIALCONDITION
#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"
#include "F77Interfaces/F77FileInitialCondition.h"

#define f_init_commongm FORTRAN_NAME(comblgm, COMBLGM)

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

#ifdef OWN_INITIALCONDITION

class InitialConditionSpecific : 
  public F77FileInitialCondition<VectorType,DIM> {
public:    
  InitialConditionSpecific() : 
    F77FileInitialCondition<VectorType,DIM>(f_initial,f_flgout,f_flgin,f_prolong) {}
};
#endif

#ifdef OWN_FLAGGING

// Threshold Flagging
template <class VectorType, class FlagType, int dim>
class F77ByThreshold : public F77OutBase<VectorType,dim>,
  public ByValue<VectorType,FlagType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  typedef ByValue<VectorType,FlagType,dim> base;
  typedef F77OutBase<VectorType,dim> out_base;
public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::grid_fct_type grid_fct_type;
  typedef typename base::flag_fct_type flag_fct_type;
  typedef typename out_base::generic_func_type generic_func_type;
  typedef typename base::max_vector_type max_vector_type;

  F77ByThreshold(generic_func_type out) : out_base(out), base() {
      for ( int i = 0 ; i < MAXCOMPONENTS ; i++ ) {
          _ValueMin[i] =  _ValueMax[i] = 0.0;
      }
  }

  virtual void register_at(ControlDevice& Ctrl) {}
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"F77ByThreshold");
    char Name[16];
    for (int d=0; d<max_vector_type::Length(); d++) {
      std::sprintf(Name,"ValueMin(%d)",d+1);
      RegisterAt(base::LocCtrl,Name,_ValueMin[d]);
      std::sprintf(Name,"ValueMax(%d)",d+1);
      RegisterAt(base::LocCtrl,Name,_ValueMax[d]);
      std::sprintf(Name,"MaxLev(%d)",d+1);
      RegisterAt(base::LocCtrl,Name,base::_MaxLevel[d]);
      std::sprintf(Name,"Type(%d)",d+1);
      RegisterAt(base::LocCtrl,Name,base::_Comparison[d]);
    }    
    RegisterAt(base::LocCtrl,"Output",base::_Output);
  }

  virtual void update() { 
    int d=max_vector_type::Length();
    
    if (d>0) {
      base::SetNcnt(d);
      base::SetIsUsed(true);
    }
    else 
      base::SetIsUsed(false);
  }

  virtual void FlagByThreshold(grid_fct_type& work, flag_fct_type& flags,
                           const int& Time, const int& Level, 
                           const DataType& ValueMin, const DataType& ValueMax, 
                           const FlagType& FlagValue,
                           const int comparison = 0, const int growby = 0) {

    int TStep = TimeStep(work,Level);
    forall(flags, Time, Level, c)
      BeginFastIndex3(work, work(Time+TStep,Level,c).bbox(), 
                      work(Time+TStep,Level,c).data(), DataType);        
      BeginFastIndex3(flags, flags(Time,Level,c).bbox(), 
                      flags(Time,Level,c).data(), FlagType);

      BBox OpBox = flags.interiorbbox(Time,Level,c);
      Coords& Os = OpBox.stepsize();
      BBox OpwBox = work(Time+TStep,Level,c).bbox();
      Coords& Ows = OpwBox.stepsize();

      for_3 (n, m, l, OpBox, Os)
        DataType val;
        val = FastIndex3(work,n-n%Ows(0),m-m%Ows(1),l-l%Ows(2));
        if ((comparison==0 && val < ValueMax && val > ValueMin) ||
            (comparison==1 && val > ValueMax && val < ValueMin)) {
          if (growby<=0) 
            FastIndex3(flags,n,m,l) = FlagValue;
          else {
            BBox bbflag(3,n,m,l,n,m,l,Os(0),Os(1),Os(2)); bbflag.grow(growby);
            flags(Time,Level,c).equals(FlagValue,bbflag);
          }
        }
      end_for

      EndFastIndex3(work);
      EndFastIndex3(flags);
    end_forall
  } 

  virtual bool SetFlags(vec_grid_fct_type& u, grid_fct_type& work, flag_fct_type& flags, 
                        const int& cnt, const int& Time, const int& Level, const double& t, 
                        const FlagType& FlagValue) {
    if ( _ValueMin[cnt] >=  _ValueMax[cnt] || 
        (base::MaxLevel(cnt)<Level && base::MaxLevel(cnt)>=0)) return false;
    int TStep = TimeStep(work,Level);
    out_base::Transform(u,work,Time,Level,cnt+1,t);
    
    forall(u,Time,Level,c)
      work(Time+TStep,Level,c).equals(work(Time,Level,c));
    end_forall

    FlagByThreshold(work, flags, Time, Level, _ValueMin[cnt], _ValueMax[cnt], 
                FlagValue, base::_Comparison[cnt]);
    return true;
  }

  virtual void OutputName(char* name, int cnt) { std::sprintf(name,"f77threshold_%d",cnt+1); }

private:
    DataType _ValueMin[MAXCOMPONENTS], _ValueMax[MAXCOMPONENTS];
};


#define NZONES 1

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));
      base::AddCriterion(new F77ByThreshold<VectorType,FlagType,DIM>(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(); }

    int LowCut(int i, int dim) { return dcmin[i][dim]; }

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();
    int dim = DIM;
    int dcmin = ((FlaggingSpecific*)_Flagging)->LowCut(0, 0);
     f_init_commongm(dim,shape,geom,dcmin);
     _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;
};