WENOIntegrator.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, Carlos Pantano, David Hill 

#ifndef AMROC_WENO_INTEGRATOR_H
#define AMROC_WENO_INTEGRATOR_H

#include "Integrator.h"
#include "CLESLog.h"

#include <string>

#ifndef f_init_common
#define f_init_common FORTRAN_NAME(combl, COMBL)
#endif

#define f_init_weno FORTRAN_NAME(initweno, INITWENO)
#define f_finish_weno FORTRAN_NAME(finishweno, FINISHWENO)

extern "C" {
  void f_init_weno(const INTEGER& dim, const INTEGER& meqn, const INTEGER& nvars, 
                   const INTEGER& mghosts, 
                   const INTEGER& order, const INTEGER& optimized,
                   const INTEGER& use_carbfix, const INTEGER& method, 
                   const INTEGER& viscous, const INTEGER& les, 
                   const INTEGER& usrc, 
                   const INTEGER& noTimeRefine, const DOUBLE& filter_alpha);
  void f_finish_weno();
  void f_init_common();
}


template <class VectorType, int dim>
class WENOIntegrator : public Integrator<VectorType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  typedef Integrator<VectorType,dim> base;

public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef generic_fortran_func generic_func_type;

  typedef int DCFlagDataType;
  typedef GridData<DCFlagDataType,dim> flag_grid_data_type;  

  typedef void (*check_1_func_type) ( FI(1,VectorType), BI, const INTEGER& meqn, 
                                      const INTEGER& mout, INTEGER& result );
  typedef void (*check_2_func_type) ( FI(2,VectorType), BI, const INTEGER& meqn, 
                                      const INTEGER& mout, INTEGER& result );
  typedef void (*check_3_func_type) ( FI(3,VectorType), BI, const INTEGER& meqn, 
                                      const INTEGER& mout, INTEGER& result );

  typedef void (*step_func_type) (const INTEGER& rk, VectorType ux[], 
                                  VectorType uxold[], const DOUBLE& dt , 
                                  DOUBLE& cfl, VectorType fxi[], 
                                  INTEGER dcflag[], INTEGER& iFilter);
  typedef void (*bnds_func_type) (INTEGER ix[], DOUBLE lbc[], DOUBLE ubc[], 
                                  DOUBLE dx[], const DOUBLE* bnd, 
                                  const INTEGER& mb, const INTEGER per[]);

  WENOIntegrator(generic_func_type step, generic_func_type bnds, generic_func_type check) : 
    base(), f_step(step), f_bnds(bnds), f_chk(check),
    _order(2), _optimized(0), _use_carbfix(1), _method(0),
    _visc(0), _les(0),_usrc(0), _check(0), _DCFlagAdaptBndry(1),
    _noTimeRefine(0), _FilterStep(-1), _FilterStrength(0.0)
  { 
    _name = ""; 
    FluxData = (VectorType*) 0;
  }

  ~WENOIntegrator() {
    if (FluxData) {
      delete [] FluxData;
      FluxData = (VectorType*) 0;
    }  
  }

  virtual double PassTimeStepFraction(int mpass) { 
    switch ( mpass ) {
    case 1:
      return 0.0;
    case 2:
      return 1.0;
    case 3:
      return 0.5;
    default:
      return 0.0;
    }
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"WENOIntegrator");
    _name = prefix+"WENOIntegrator";
    RegisterAt(LocCtrl,"Order",_order);
    RegisterAt(LocCtrl,"Optimized",_optimized);
    RegisterAt(LocCtrl,"CarbuncleFix",_use_carbfix);
    RegisterAt(LocCtrl,"Method",_method);
    RegisterAt(LocCtrl,"UseViscous",_visc);
    RegisterAt(LocCtrl,"UseLES",_les);
    RegisterAt(LocCtrl,"UseSource",_usrc);
    RegisterAt(LocCtrl,"Check",_check);
    RegisterAt(LocCtrl,"DCFlagAdaptBndry",_DCFlagAdaptBndry);
    
    RegisterAt(LocCtrl,"FilterStep",_FilterStep);
    RegisterAt(LocCtrl,"FilterStrength",_FilterStrength);

    log.register_at(LocCtrl, ""); 
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }
   
  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
    base::SetupData(gh,ghosts);
    base::SetMaxIntegratorPasses(NMaxPass());
    base::_abort = (NCheck()>=10 ? DAGHTrue : DAGHFalse);
    int _dim = dim;

    log.initialize();

    f_init_weno(_dim, base::NEquations(), NVARS, base::NGhosts(),
                _order, _optimized, _use_carbfix, _method, _visc, _les,
                _usrc, _noTimeRefine, _FilterStrength);
    f_init_common();
  }

  virtual void finish() { f_finish_weno(); log.finalize(); } 

  virtual double CalculateGrid(vec_grid_data_type& NewStateVec, 
                               vec_grid_data_type& OldStateVec,
                               vec_grid_data_type* Flux[],
                               const int& level, const double& t, 
                               const double& dt, const int& mpass) {
    double cfl = 0.0;
#ifdef DEBUG_PRINT_INTEGRATOR
    ( comm_service::log() << "on: " << OldStateVec.bbox() ).flush();
#endif

    if (mpass==1 || mpass==0) NewStateVec.copy(OldStateVec);

    int mx[dim], d, per[dim];
    Coords ex = NewStateVec.extents();
    for (d=0; d<dim; d++) { 
      mx[d] = ex(d)-2*base::NGhosts();
      per[d] = base::GH().periodicboundary(d);
    }

    flag_grid_data_type** DCFlag;
    DCFlagDataType* DCFlagData;
    AllocDCFlag(NewStateVec.bbox(), DCFlag, DCFlagData);
    if (mpass==1 || mpass==0) SetDCFlag(NewStateVec.bbox(), DCFlag, level);

    DCoords dx = base::GH().worldStep(NewStateVec.stepsize());
    DCoords lbc = base::GH().worldCoords(NewStateVec.lower(), NewStateVec.stepsize());
    DCoords ubc = base::GH().worldCoords(NewStateVec.upper(), NewStateVec.stepsize());

    ((bnds_func_type) f_bnds)(mx,lbc(),ubc(),dx(),base::GH().wholebndry(),
                              base::GH().nbndry(),per);
    
    if (NCheck()>0 && (mpass==1 || mpass==0))
      base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::before f_step");

    int iFilter = 0;
    int Time = CurrentTime(base::GH(), 0 );
    if ( _FilterStep>0 && Time%(_FilterStep*StepSize(base::GH(),0)) == 0 ) iFilter = 1;

    // The ghostcells must remain untouched during integration!
    ((step_func_type) f_step)(mpass,NewStateVec.data(),OldStateVec.data(), 
                              dt,cfl,Flux[0]->data(),DCFlag[0]->data(),iFilter);

    if (NCheck()>=0) 
      // The negation captures cfl=NaN!
      if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0 || !(cfl>=0) ||
          (!(cfl>0) && (mpass==NMaxPass()))) {  
        ResetGridFluxes(Flux);
        NewStateVec.copy(OldStateVec);
        if (cfl>0) cfl = std::max(cfl, 2.0);
      }
    
    DeAllocDCFlag(DCFlag,DCFlagData);

#ifdef DEBUG_PRINT_INTEGRATOR
      ( comm_service::log() << "  CFL = " << cfl <<
        "   t = " << t << "   dt = " << dt << "\n" ).flush();
#endif
    
    return cfl;
  } 

  void SetDCFlag(const BBox &databbox, flag_grid_data_type** DCFlag, const int& l) {
    int d;
    for (d=0; d<dim; d++) DCFlag[d]->equals(0);

    if ( _DCFlagAdaptBndry==0 ) return;

    BBox interiorbbox = grow(databbox, -base::NGhosts());

    // check if need to flag fine side.
    if (base::U().prolongbndrylevel(l)!=0 ) {
      for (int s=0; s<2*dim; s++) { 
        d = s/2;
        BBox localbbndry(interiorbbox);
        if (s%2 == 0) {
          localbbndry.growlower(d,-1); 
          localbbndry.upper(d) = localbbndry.lower(d);
        }
        else {
          localbbndry.growupper(d,1); 
          localbbndry.lower(d) = localbbndry.upper(d);
        }
        BBoxList localbndrylist = *base::U().prolongbndrylevel(l);
        localbndrylist *= localbbndry;
        for (register BBox* bb=localbndrylist.first();bb;
             bb=localbndrylist.next()) {
          if (s%2 == 0) {
            bb->lower(d) += bb->stepsize(d);
            bb->upper(d) = bb->lower(d);
          } 
          DCFlag[d]->equals((-1+2*(s%2)),*bb);
        }
      }
    }

    if ( _DCFlagAdaptBndry==1 ) return;
    
    // check if we need to flag coarse side
    if ( l < FineLevel(base::GH()) ) {
      int by = StepSize(base::GH(),l)/StepSize(base::GH(),l+1);
#ifdef _DEBUG_DCFLAG
      printf("\nGhost cells %d\n", base::NGhosts());
      printf("father(%d) (%d %d) (%d %d) [%d %d]\n",l,
             interiorbbox.lower(0), interiorbbox.upper(0), 
             interiorbbox.lower(1), interiorbbox.upper(1), 
             interiorbbox.stepsize(0), interiorbbox.stepsize(1));
#endif
      GridBoxList* gbchild = base::GH().ggbl(l+1);
      if ( gbchild )
        for (register GridBox *g=gbchild->first();g;g=gbchild->next())
          {
            //BBox bb = grow(g->gbBBox(), base::NGhosts());
            BBox bb(g->gbBBox());
            bb.coarsen(by);
#ifdef _DEBUG_DCFLAG
            printf("child(%d) (%d %d) (%d %d) [%d %d]\n",l+1,
                   bb.lower(0), bb.upper(0), 
                   bb.lower(1), bb.upper(1), 
                   bb.stepsize(0), bb.stepsize(1));
#endif
            
            for (int s=0; s<2*dim; s++) { 
              d = s/2;
              BBox localbbndry(bb);
              if (s%2 == 0) {
                localbbndry.growlower(d,-1); 
                localbbndry.upper(d) = localbbndry.lower(d);
              }
              else {
                localbbndry.growupper(d,1); 
                localbbndry.lower(d) = localbbndry.upper(d);
              }
#ifdef _DEBUG_DCFLAG
              printf("to intersect[%d] (%d %d) (%d %d) [%d %d]\n", s,
                     localbbndry.lower(0), localbbndry.upper(0), 
                     localbbndry.lower(1), localbbndry.upper(1), 
                     localbbndry.stepsize(0), localbbndry.stepsize(1));
#endif
              localbbndry *= interiorbbox;
              if ( ! localbbndry.empty() ) {
                if ( s%2 == 0 ) {
                  localbbndry.growlower(d,1); 
                  localbbndry.upper(d) = localbbndry.lower(d);
                }
#ifdef _DEBUG_DCFLAG
                printf("intersection (%d %d) (%d %d) [%d %d]\n",
                       localbbndry.lower(0), localbbndry.upper(0), 
                       localbbndry.lower(1), localbbndry.upper(1), 
                       localbbndry.stepsize(0), localbbndry.stepsize(1));
#endif
                
                DCFlag[d]->equals((1-2*(s%2)),localbbndry);
              }
            }
          }
    }
  }
    
  void AllocDCFlag(const BBox &databbox, flag_grid_data_type**& DCFlag, 
                   DCFlagDataType*& DCFlagData) {
    BBox fluxbbox = grow(databbox, -base::NGhosts());
    fluxbbox.growupper(1);
    // Create a CONSECUTIVE array for the flags in ALL directions!
    DCFlagData = new DCFlagDataType[base::Dim()*fluxbbox.size()];
    DCFlag = new flag_grid_data_type* [base::Dim()]; 
    for (register int d=0; d<base::Dim(); d++)
      DCFlag[d] = new flag_grid_data_type(fluxbbox,&(DCFlagData[d*fluxbbox.size()])); 
  }

  void DeAllocDCFlag(flag_grid_data_type**& DCFlag, DCFlagDataType*& DCFlagData) {
    if (DCFlag) {
      DCFlagDataType* dummy;
      for (register int d=0; d<base::Dim(); d++) 
        if (DCFlag[d]) {
          // Erease pointer to data
          DCFlag[d]->deallocate(dummy);
          delete DCFlag[d];
        }
      delete [] DCFlag;
    }
    if (DCFlagData) {
      delete [] DCFlagData;
      DCFlagData = (DCFlagDataType*) 0;
    }
  }
    
  virtual void AllocGridFluxes(const BBox &bb, vec_grid_data_type**& Flux) {
    // The solver assumes a CONSECUTIVE array for the fluxes in ALL directions!
    // We need a member variable to store it between different member function calls.
    FluxData = new VectorType[base::Dim()*bb.size()];
    Flux = new vec_grid_data_type* [2*base::Dim()]; 
    for (register int d=0; d<base::Dim(); d++) {
      Flux[2*d] = new vec_grid_data_type(bb,&(FluxData[d*bb.size()])); 
      Flux[2*d+1] = Flux[2*d];
    }
  }

  virtual void DeAllocGridFluxes(vec_grid_data_type**& Flux) {
    if (Flux) {
      VectorType* dummy;
      for (register int d=0; d<base::Dim(); d++) 
        if (Flux[2*d]) {
          // Erease pointer to data
          Flux[2*d]->deallocate(dummy);
          delete Flux[2*d];
        }
      delete [] Flux;
    }
    if (FluxData) {
      delete [] FluxData;
      FluxData = (VectorType*) 0;
    }
  }
    
  virtual void ResetGridFluxes(vec_grid_data_type**& Flux) {
    if (Flux) {
      VectorType zero(0.0);
      for (register int d=0; d<base::Dim(); d++) 
        if (Flux[2*d]) 
          Flux[2*d]->equals(zero);
    }
  }

  virtual int ControlGrid(vec_grid_data_type& StateVec, const int& level, 
                          const BBox& where, const double& time, const int verbose) {
    int result = 1;
    if (!f_chk) return result;
    if (dim == 1) 
      ((check_1_func_type) f_chk)(FA(1,StateVec),BOUNDING_BOX(where),base::NEquations(),
                                  verbose,result);
    else if (dim == 2) 
      ((check_2_func_type) f_chk)(FA(2,StateVec),BOUNDING_BOX(where),base::NEquations(),
                                  verbose,result);
    else if (dim == 3) 
      ((check_3_func_type) f_chk)(FA(3,StateVec),BOUNDING_BOX(where),base::NEquations(),
                                  verbose,result);
    return result; 
  }  

  inline const int& NCheck() const { return _check; }
  virtual int NMethodOrder() const { return _order; }
  inline int NMaxPass() const { return 3; }
  virtual void SetNoTimeRefine(int flag) { _noTimeRefine = flag; }

  inline void SetCheckFunc(generic_func_type check) { f_chk = check; }
  generic_func_type GetCheckFunc() const { return f_chk; }
  inline void SetStepFunc(generic_func_type step) { f_step = step; }
  generic_func_type GetStepFunc() const { return f_step; }
  inline void SetBndsFunc(generic_func_type bnds) { f_bnds = bnds; }
  generic_func_type GetBndsFunc() const { return f_bnds; }

protected:
  generic_func_type f_step, f_bnds, f_chk;
  std::string _name;
  int _order,_optimized, _use_carbfix;
  int _method, _visc, _les, _usrc, _check;
  int _DCFlagAdaptBndry;
  int _noTimeRefine;
  VectorType* FluxData;
  CLESLog log;
  int _FilterStep;
  double _FilterStrength;
};



#endif