ClpIntegratorBase.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus
//
// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_CLP_INTEGRATORBASE_H
#define AMROC_CLP_INTEGRATORBASE_H

#include "Integrator.h"
#include "Timing.h"
#include "FixupBase.h"

#include <iostream>
#include <string>
#include <cfloat>

#define FACTOR 1.0e5

#define f_init_common FORTRAN_NAME(combl, COMBL)
extern "C" {
   void f_init_common();
}

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

public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef GridData<VectorType,minus_1<dim>::dim>    ld_vec_grid_data_type;
  typedef GridData<AuxVectorType,minus_1<dim>::dim> ld_aux_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 (*aux_1_func_type) ( const INTEGER& maxmx,
                                    const INTEGER& meqn, const INTEGER& mbc,
                                    const INTEGER& ibx,
                                    const INTEGER& mx, VectorType q[],
                                    DOUBLE aux[], const INTEGER& maux,
                                    const DOUBLE& cornx,
                                    const DOUBLE& dx,
                                    const DOUBLE& t, const DOUBLE& dt );
  typedef void (*aux_2_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy,
                                    const INTEGER& meqn, const INTEGER& mbc,
                                    const INTEGER& ibx, const INTEGER& iby,
                                    const INTEGER& mx, const INTEGER& my, VectorType q[],
                                    DOUBLE aux[], const INTEGER& maux,
                                    const DOUBLE& cornx, const DOUBLE& corny,
                                    const DOUBLE& dx, const DOUBLE& dy,
                                    const DOUBLE& t, const DOUBLE& dt );
  typedef void (*aux_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
                                    const INTEGER& meqn, const INTEGER& mbc,
                                    const INTEGER& ibx, const INTEGER& iby, const INTEGER& ibz,
                                    const INTEGER& mx, const INTEGER& my, const INTEGER& mz, VectorType q[],
                                    DOUBLE aux[], const INTEGER& maux,
                                    const DOUBLE& cornx, const DOUBLE& corny, const DOUBLE& cornz,
                                    const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz,
                                    const DOUBLE& t, const DOUBLE& dt );

  typedef void (*source_1_func_type) ( const INTEGER& maxmx,
                                       const INTEGER& meqn, const INTEGER& mbc,
                                       const INTEGER& ibx,
                                       const INTEGER& mx, 
                                       VectorType q[],
                                       const DOUBLE aux[], const INTEGER& maux,
                                       const DOUBLE& t, const DOUBLE& dt, const INTEGER& ibnd );
  typedef void (*source_2_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy,
                                       const INTEGER& meqn, const INTEGER& mbc,
                                       const INTEGER& ibx, const INTEGER& iby,
                                       const INTEGER& mx, const INTEGER& my, 
                                       VectorType q[],
                                       const DOUBLE aux[], const INTEGER& maux,
                                       const DOUBLE& t, const DOUBLE& dt, const INTEGER& ibnd );
  typedef void (*source_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy,  const INTEGER& maxmz,
                                       const INTEGER& meqn, const INTEGER& mbc,
                                       const INTEGER& ibx, const INTEGER& iby, const INTEGER& ibz,
                                       const INTEGER& mx, const INTEGER& my, const INTEGER& mz, 
                                       VectorType q[],
                                       const DOUBLE aux[], const INTEGER& maux,
                                       const DOUBLE& t, const DOUBLE& dt, const INTEGER& ibnd );

  ClpIntegratorBase(const int equ, const int wv) : 
    base(), _EqUsed(equ), _Waves(wv), f_chk(0), f_aux(0), f_src(0)
  { ConsInit(); }
  ClpIntegratorBase(const int equ, const int wv, generic_func_type check, 
                    generic_func_type aux, generic_func_type source) : 
    base(), _EqUsed(equ), _Waves(wv), f_chk(check), f_aux(aux), f_src(source)
  { ConsInit(); }

  void ConsInit() {
    limiter = new int[NWaves()];
    mthlim = new int[NWaves()];
    for (int i=0; i<NWaves(); i++) {
      limiter[i] = -1; mthlim[i] = -1;
    }
    _name = "";
    for (int d=0; d<7; d++) method[d] = 0;
  }    

  ~ClpIntegratorBase() {
    if (limiter) delete [] limiter;
    if (mthlim) delete [] mthlim;
  }

  //******************************************************************************
  // Abstract class interface
  //******************************************************************************
  virtual void ComputeRP(ld_vec_grid_data_type& LeftState, ld_aux_grid_data_type& auxl, 
                         ld_vec_grid_data_type& RightState, ld_aux_grid_data_type& auxr, 
                         ld_vec_grid_data_type& NewState, const int& direction) = 0;
  virtual double ComputeGrid(vec_grid_data_type& StateVec, 
                             vec_grid_data_type* Flux[], DataType* aux,
                             const double& t, const double& dt, const int& mpass) = 0;
  //******************************************************************************

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"ClawpackIntegrator");
    _name = prefix+"ClawpackIntegrator";
    char MethodName[16];
    register int d;
    for (d=0; d<7; d++) {
      sprintf(MethodName,"Method(%d)",d+1);
      RegisterAt(LocCtrl,MethodName,method[d]);
    }
    char LimiterName[16];
    for (d=0; d<NWaves(); d++) {
      sprintf(LimiterName,"Limiter(%d)",d+1);
      RegisterAt(LocCtrl,LimiterName,limiter[d]);
    }
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }
   
  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
    for (int i=0; i<NWaves(); i++) 
      if (limiter[i] > 0) mthlim[i] = limiter[i];
      else mthlim[i] = limiter[0];
    if (NAux()!=0 && NAux()!=AuxVectorType::Length())
      method[6] = 0;
    base::SetMaxIntegratorPasses(NMaxPass());
    base::_abort = (NCheck()>=10 ? DAGHTrue : DAGHFalse);
    f_init_common();
    base::SetupData(gh,ghosts);
  }

  virtual void finish() {
    if (limiter) {
      delete [] limiter;
      limiter = (int *) 0;
    }
    if (mthlim) {
      delete [] mthlim;
      mthlim = (int *) 0;
    }
    base::finish();
  } 

  void CalculateRP(ld_vec_grid_data_type &LeftState, BBox& bbl, double tl, double dtl, 
                   ld_vec_grid_data_type &RightState, BBox& bbr, double tr, double dtr, 
                   ld_vec_grid_data_type &NewState, const int& direction) {
    
    //------------------------------------------------------------------                      
    // Calculate AuxArray and source term on left and right state.
    //------------------------------------------------------------------                      
    vec_grid_data_type left(bbl, LeftState.data());
    vec_grid_data_type right(bbr, RightState.data());
    AllocError::SetTexts("ClpIntegrator","calculate_Riemann_Problem(): aux arrays");
    DataType* auxldata = new DataType[LeftState.bbox().size()*NAux()];
    DataType* auxrdata = new DataType[RightState.bbox().size()*NAux()];

    // Set aux arrays.
    if (NAux()>0) {
      SetAux(auxldata,left,tl,dtl);
      SetAux(auxrdata,right,tr,dtr);
    }      

    // Apply the source term, if necessary.
    START_INTERMEDIATE_WATCH
    START_WATCH
      if (NSrc() > 0) {
        // Integrate source term on left and right data to the same time level.
        int srcbnd = 1;
        if (tl < tr) {
          SetSrc(auxldata,left,tl,tr-tl,srcbnd);
          tl += tr-tl;
        }
        if (tl > tr) {
          SetSrc(auxrdata,right,tr,tl-tr,srcbnd);
          tr += tl-tr;
        }

        if (NSrc() > 1) {
          double dt = Min(dtl,dtr);
          if (NSrc()==2) {
            SetSrc(auxldata,left,tl,dt/2.0,srcbnd);
            SetSrc(auxrdata,right,tr,dt/2.0,srcbnd);
          }
          if (NSrc()==3) {
            SetSrc(auxldata,left,tl,dt,srcbnd);
            SetSrc(auxrdata,right,tr,dt,srcbnd);
          }
          if (NAux()>0) {
            SetAux(auxldata,left,tl,dtl);
            SetAux(auxrdata,right,tr,dtr);
          }      
        }         
      }
    END_WATCH_SOURCE_INTEGRATION
    END_INTERMEDIATE_WATCH
    VectorType* data_dummy;
    left.deallocate(data_dummy);
    right.deallocate(data_dummy);
    ld_aux_grid_data_type auxl(LeftState.bbox(), (AuxVectorType*) auxldata);
    ld_aux_grid_data_type auxr(RightState.bbox(), (AuxVectorType*) auxrdata);

    // Solve the Riemann problems. 
    ComputeRP(LeftState, auxl, RightState, auxr, NewState, direction);

    AuxVectorType* aux_data_dummy;
    auxl.deallocate(aux_data_dummy);
    auxr.deallocate(aux_data_dummy);
    delete [] auxldata; delete [] auxrdata;
  }

  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 (NMaxPass()==0 || mpass==1) 
      NewStateVec.copy(OldStateVec);

    if (dt >= DBL_EPSILON*FACTOR) {
      AllocError::SetTexts("ClpIntegrator","calculate_time_step(): aux");
      DataType* aux = new DataType[OldStateVec.bbox().size()*NAux()];     
      if (NAux() > 0) SetAux(aux,NewStateVec,t,dt);

#ifdef DEBUG_AMR
      if (NCheck()%10>1 && (NSrc()==2 || NSrc()==3))
        base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::before SrcInt");
#endif

      // The source term has to be applied also on ghostcells, otherwise boundary conditions
      // would not be consistent. See CLAWPACK-Manual Note #5 "Source Terms".
      START_INTERMEDIATE_WATCH
      START_WATCH
        if (NSrc()>1 && (NMaxPass()==0 || mpass==1)) {
          int srcbnd = 1;
          if (NSrc()==2) SetSrc(aux,NewStateVec,t,dt/2.0,srcbnd);
          if (NSrc()==3) SetSrc(aux,NewStateVec,t,dt,srcbnd);
          if (NAux()>0) SetAux(aux,NewStateVec,t,dt);
          if (NCheck()>=0) if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0) {
            ResetGrid(NewStateVec,OldStateVec,Flux,cfl,level,t,mpass);
            delete [] aux;
            if (cfl>0) cfl = std::max(cfl, 2.0);
            return cfl;
          }
        }
      END_WATCH_SOURCE_INTEGRATION
      END_INTERMEDIATE_WATCH

      if (NCheck()>0 && (NMaxPass()==0 || mpass==1))
        base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::before f_step");

      // Update NewStateVec. 
      cfl = ComputeGrid(NewStateVec, Flux, aux, t, dt, mpass);

      if (NCheck()>=0) 
        // The negation captures cfl=NaN!
        if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0 || !(cfl>0)) {
          ResetGrid(NewStateVec,OldStateVec,Flux,cfl,level,t,mpass);
          delete [] aux;
          if (cfl>0) cfl = std::max(cfl, 2.0);
          return cfl;
        }

#ifdef DEBUG_AMR
      if (NCheck()%10>1 && (NMaxPass()==0 || mpass==NMaxPass()))
        base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::after f_step");
#endif

      // The source term is not applied to ghostcells, because they are overwritten anyway by
      // the application of boundary conditions in the next step.
      START_INTERMEDIATE_WATCH
      START_WATCH
        if (NSrc()>0 &&  NSrc()<3 && (NMaxPass()==0 || mpass==NMaxPass())) {
          int srcbnd = 0;
          if (NAux()>0) SetAux(aux,NewStateVec,t,dt);
          if (NSrc()==2) SetSrc(aux, NewStateVec, t+dt/2.0, dt/2.0, srcbnd);
          if (NSrc()==1) SetSrc(aux, NewStateVec, t, dt, srcbnd);
          if (NCheck()>=0) if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0) {
            ResetGrid(NewStateVec,OldStateVec,Flux,cfl,level,t,mpass);
            delete [] aux;
            if (cfl>0) cfl = std::max(cfl, 2.0);
            return cfl;
          }
        }
      END_WATCH_SOURCE_INTEGRATION
      END_INTERMEDIATE_WATCH

#ifdef DEBUG_AMR
      if (NCheck()%10>1 && (NSrc()==2 || NSrc()==1) && (mpass==0 || mpass==NMaxPass()))
        base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::after SrcInt");
#endif

      delete [] aux;
    }
    else {
      VectorType zero(0.0);
      for (int d=0; d<dim; d++)
        Flux[d]->equals(zero);
#ifdef DEBUG_PRINT_INTEGRATOR
      ( comm_service::log() <<  "dt < eps*1.0e5  " ).flush();
#endif
    }

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

  virtual void AllocGridFluxes(const BBox &bb, vec_grid_data_type**& Flux) {
    Flux = new vec_grid_data_type* [2*base::Dim()]; 
    for (register int d=0; d<2*base::Dim(); d++) 
      Flux[d] = new vec_grid_data_type(bb); 
  }

  virtual void DeAllocGridFluxes(vec_grid_data_type**& Flux) {
    if (Flux) {
      for (register int d=0; d<2*base::Dim(); d++) 
        if (Flux[d]) delete Flux[d];
      delete [] Flux;
    }
  }
    
  virtual void ResetGridFluxes(vec_grid_data_type**& Flux) {
    if (Flux) {
      VectorType zero(0.0);
      for (register int d=0; d<2*base::Dim(); d++) 
        if (Flux[d]) 
          Flux[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; 
  }  

  void ResetGrid(vec_grid_data_type& StateVec, vec_grid_data_type& RecoverStateVec,
                 vec_grid_data_type* Flux[], double& cfl, const int& level, 
                 const double& t, const int& mpass) {
    ResetGridFluxes(Flux);
    StateVec.copy(RecoverStateVec);
    if (cfl>0) cfl = std::max(cfl, 2.0);
    //#ifdef DEBUG_PRINT_INTEGRATOR
#ifdef DEBUG_PRINT
    ( comm_service::log() << "Data in " << StateVec.bbox() 
      << " recovered (Level " << level << "  Time = " << t << "  mpass = " 
      << mpass << "). Using CFL=" << cfl << "." << std::endl ).flush();
#endif
  }

  void SetAux(DataType* aux, vec_grid_data_type& StateVec, const double& t, const double& dt) {
    if (NAux() <=0 || !f_aux) return;    
    DCoords dx = base::GH().worldStep(StateVec.stepsize());
    DCoords lbcorner = base::GH().worldCoords(StateVec.lower(),StateVec.stepsize());
    int mx[3], ibx[3];
    DimExtBBox(StateVec.bbox(), mx, ibx);
    if (dim == 1)
      ((aux_1_func_type) f_aux)(AA(1,mx),base::NEquations(),base::NGhosts(),AA(1,ibx),AA(1,mx),
                                StateVec.data(),aux,NAux(),AA(1,lbcorner()),AA(1,dx()),t,dt);   
    else if (dim == 2)
      ((aux_2_func_type) f_aux)(AA(2,mx),base::NEquations(),base::NGhosts(),AA(2,ibx),AA(2,mx),
                                StateVec.data(),aux,NAux(),AA(2,lbcorner()),AA(2,dx()),t,dt); 
    else if (dim == 3)
      ((aux_3_func_type) f_aux)(AA(3,mx),base::NEquations(),base::NGhosts(),AA(3,ibx),AA(3,mx),
                                StateVec.data(),aux,NAux(),AA(3,lbcorner()),AA(3,dx()),t,dt); 
  }

  virtual void SetSrc(DataType* aux, vec_grid_data_type& StateVec, const double t, 
                       const double dt, const int& bnd) {
    if (NSrc() <=0 || !f_src) return;    
    DCoords dx = base::GH().worldStep(StateVec.stepsize());
    DCoords lbcorner = base::GH().worldCoords(StateVec.lower(),StateVec.stepsize());
    int mx[3], ibx[3];
    DimExtBBox(StateVec.bbox(), mx, ibx);
    if (dim == 1)
      ((source_1_func_type) f_src)(AA(1,mx),base::NEquations(),base::NGhosts(),AA(1,ibx),AA(1,mx),
                                   StateVec.data(),aux,NAux(),t,dt,bnd); 
    else if (dim == 2)
      ((source_2_func_type) f_src)(AA(2,mx),base::NEquations(),base::NGhosts(),AA(2,ibx),AA(2,mx),
                                   StateVec.data(),aux,NAux(),t,dt,bnd); 
    else if (dim == 3)
      ((source_3_func_type) f_src)(AA(3,mx),base::NEquations(),base::NGhosts(),AA(3,ibx),AA(3,mx),
                                   StateVec.data(),aux,NAux(),t,dt,bnd); 
  }

  void DimExtBBox(const BBox& bb, int mx[], int ibx[]) {
    for (int d=0; d<dim; d++) 
      if (bb.extents(d)-2*base::NGhosts()<1) 
        { mx[d] = bb.extents(d); ibx[d] = 0; }
      else
        { mx[d] = bb.extents(d)-2*base::NGhosts(); ibx[d] = 1; }
  }

  virtual int NMethodOrder() const { return (std::max(method[1],2)); }
  inline const int& NCheck() const { return method[3]; }
  inline int NSrc() const { return (method[4]%10); }
  inline bool DimSplit() const { return (method[2]<0); }
  inline int NMaxPass() const { return (method[2]<-1 && dim>1 ? dim : 0); }
  inline const int& NAux() const { return method[6]; }
  inline void SetEqUsed(const int& eq) { _EqUsed = eq; }
  inline const int& NEqUsed() const { return _EqUsed; }
  inline void SetWaves(const int& wv) { _Waves = wv; }
  inline const int& NWaves() const { return _Waves; }  

  inline void SetCheckFunc(generic_func_type check) { f_chk = check; }
  generic_func_type GetCheckFunc() const { return f_chk; }
  inline void SetAuxFunc(generic_func_type aux) { f_aux = aux; }
  generic_func_type GetAuxFunc() const { return f_aux; }
  inline void SetSrcFunc(generic_func_type src) { f_src = src; }
  generic_func_type GetSrcFunc() const { return f_src; }

  std::ostream& debug(std::ostream& out, const self& P) {
    out << "ClpIntegrator registered at " << P._name << "\n";
    out << "   Methods:";
    for (int i=0; i<7; i++)
      out << " " << P.method[i];
    out << "   Limiter:";
    for (int j=0; j<P.NWaves(); j++)
      out << " " << P.limiter[j];
    out << "\n";
    return out;
  }

protected:
  int _EqUsed, _Waves;
  int method[7];
  int *limiter, *mthlim;
  std::string _name;
  generic_func_type f_chk, f_aux, f_src;
  int order;
};


template <class VectorType, class AuxVectorType, int dim> class ClpIntegrator;

#endif