RIMIntegrator.h

Go to the documentation of this file.

// -*- C++ -*-
// $Header: /home/proj/vtf3d/vtf/amroc/rim/RIMIntegrator.h,v 1.6 2006/05/27 05:59:45 ralf Exp $

#ifndef AMROC_RIM_INTEGRATOR_H
#define AMROC_RIM_INTEGRATOR_H

#include "Integrator.h"

#include <string>

#define f_init_common FORTRAN_NAME(combl, COMBL)

extern "C" {
  void f_init_common();
}

template <class VectorType, int dim>
class RIMIntegrator : 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 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& nx, const INTEGER& ny, const INTEGER& nbc, 
                                  const INTEGER& neqn, const DOUBLE& dx,
                                  const DOUBLE& dy, VectorType ux[], VectorType uxold[], 
                                  const DOUBLE& dt , DOUBLE& cfl, VectorType fx[], VectorType fy[]);

  RIMIntegrator(generic_func_type step) : 
    base(), f_step(step), f_chk(0),
    _order(2), _check(0) { 
    _name = ""; 
    FluxData = (VectorType*) 0;
  }

  RIMIntegrator(generic_func_type step, generic_func_type chk) : 
    base(), f_step(step), f_chk(chk),
    _order(2), _check(0) { 
    _name = ""; 
    FluxData = (VectorType*) 0;
  }

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

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"RIMIntegrator");
    _name = prefix+"RIMIntegrator";
    RegisterAt(LocCtrl,"Check",_check);
  }
  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);
    f_init_common();
  }

  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;

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

    if (NCheck()>0)
      base::CheckGrid(OldStateVec, level, OldStateVec.bbox(), t, "CalculateGrid::before f_step");

    // The ghostcells must remain untouched during integration!
    ((step_func_type) f_step)(mx[0],mx[1],base::NGhosts(),base::NEquations(),
                              dx(0),dx(1),NewStateVec.data(),OldStateVec.data(), 
                              dt,cfl,Flux[0]->data(),Flux[2]->data());
    
    if (NCheck()>=0) 
      // The negation captures cfl=NaN!
      if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0 || !(cfl>0)) {
        ResetGridFluxes(Flux);
        NewStateVec.copy(OldStateVec);
        if (cfl>0) cfl = std::max(cfl, 2.0);
        return cfl;
      }

    if (NCheck()%10>1)
      base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::after f_step");

    return cfl;
  } 

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

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

protected:
  generic_func_type f_step, f_chk;
  std::string _name;
  int _order, _check;
  VectorType* FluxData;
};

#endif