amroc/lbm/LBMFixup.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2013 Ralf Deiterding
// German Aerospace Center (DLR)

#ifndef AMROC_LBM_FIXUP_H
#define AMROC_LBM_FIXUP_H

#include "AMRFixup.h"
#include "LBMFixupOps.h"

template <class LBMType, class FixupType, int dim>
class LBMFixup : public AMRFixup<typename LBMType::MicroType,FixupType,dim>,
                 public LBMFixupOps<LBMType,dim> {
  typedef AMRFixup<typename LBMType::MicroType,FixupType,dim> base;
  typedef AMRBase<typename LBMType::MicroType,dim> amrbase;
  typedef LBMFixupOps<LBMType,dim> lbmfixupops;

public:
  typedef typename LBMType::MicroType MicroType;
  typedef typename LBMType::MacroType MacroType;
  typedef typename LBMType::DataType DataType;

  typedef typename base::vec_grid_data_type vec_grid_data_type;
  
  LBMFixup(LBMType &lbm) : _lbm(lbm), _RescaleNonEq(0) {
    for (int d=0; d<2*dim; d++) 
      indices[d] = (int *) 0;
  }
  virtual ~LBMFixup() {
    for (int d=0; d<2*dim; d++) 
      if (indices[d]) delete [] indices[d];
  } 
    
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"Fixup");
    RegisterAt(base::LocCtrl,"RescaleNonEq",_RescaleNonEq);
  }
  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl,""); }

  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
    amrbase::SetupData(gh, ghosts);
    for (int d=0; d<2*dim; d++) {
      indices[d] = new int[LBM().NMicroVar()];
      Nindices[d] = LBM().OutgoingIndices(d,indices[d]);
      for (register int n=Nindices[d]; n<LBM().NMicroVar(); n++) 
        indices[d][n]=-1;
    }
  }

  virtual void SaveFluxes(const int Time, const int Level, const int c,
                          vec_grid_data_type* flux[], 
                          const double dt, const int& mpass) {}

  virtual void AddIntercellFluxes(const int Time, const int Level, const int c,
                                  vec_grid_data_type* flux[], const double dt,
                                  const int& mpass) {}
  
  virtual void Correction(const int Time, const int WTime, const int Level, 
                          const double t, const double dt) {

    int TimeF = CurrentTime(base::GH(),Level+1);
    int TimeC = CurrentTime(base::GH(),Level);
    const int refinefactor = RefineFactor(base::GH(),Level);
    DataType fac = refinefactor*LBM().Omega(dt/refinefactor)/LBM().Omega(dt);

    forall (base::U(),TimeF,Level+1,c) 
      for (int d=0; d<2*dim; d++) {
        if (base::ValidSide(TimeF,Level+1,c,d)) {
          BBox bb = base::U()(TimeF,Level+1,c).bbox();
          BBox from(CorrectionBBox(coarsen(base::U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
          if (d%2 == 0) from.growupper(d/2,1);
          else from.growlower(d/2,-1);                    
          BBox to(base::FixupBBox(coarsen(base::U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
          vec_grid_data_type uh(from); 
          lbmfixupops::average_indices(uh,base::U()(TimeF,Level+1,c),from,indices[d],Nindices[d]); 
          LBM().ReverseStream(uh,to,d);
          for (register int j=0; j<base::U().parents(TimeF,Level+1,c); j++) {
            BBox wto(coarsen(base::U().parentbox(TimeF,Level+1,c,j),refinefactor)*to);
            if (!wto.empty()) {
              const int& idx = base::U().parentidx(TimeF,Level+1,c,j);
              if (_RescaleNonEq) 
                lbmfixupops::scale_indices(base::U()(TimeC,Level,idx),uh,wto,indices[d],Nindices[d],fac); 
              else
                lbmfixupops::copy_indices(base::U()(TimeC,Level,idx),uh,wto,indices[d],Nindices[d]);            
            }
          }
        }
      }
    end_forall 

    END_WATCH_FIXUP_WHOLE
    START_WATCH
      AdaptiveBoundaryUpdate(base::U(),TimeC,Level);
    END_WATCH(BOUNDARIES_INTERPOLATION) 
    START_WATCH
      Sync(base::U(),TimeC,Level);
    END_WATCH_FIXUP_SYNC
    START_WATCH
      ExternalBoundaryUpdate(base::U(),TimeC,Level);
    END_WATCH(BOUNDARIES_EXTERNAL)
    START_WATCH

    forall (base::U(),TimeF,Level+1,c) 
      for (int d=0; d<2*dim; d++) {
        if (base::ValidSide(TimeF,Level+1,c,d)) {
          BBox corbb(CorrectionBBox(coarsen(base::U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
          for (register int j=0; j<base::U().parents(TimeF,Level+1,c); j++) {
            const int& idx = base::U().parentidx(TimeF,Level+1,c,j);
            BBox to = base::U().interiorbbox(Time,Level,idx)*corbb;
            if (!to.empty())
              LBM().LocalStep(base::U()(Time,Level,idx),base::U()(TimeC,Level,idx),to,dt);
          }
        }
      }
    end_forall 
  }

  BBox CorrectionBBox(const BBox& interior, const int s) {    
    int d = s/2;
#ifdef DEBUG_PRINT
    assert (d>=0 && d<dim);
#endif
    BBox bb(interior); 
    bb.grow(1);
    if (s%2 == 0) 
      bb.upper(d) = bb.lower(d);
    else 
      bb.lower(d) = bb.upper(d);
    return bb;
  }

  inline LBMType& LBM() { return _lbm; }
  inline const LBMType& LBM() const { return _lbm; }

protected:
  LBMType &_lbm;
  int _RescaleNonEq;
  int Nindices[2*dim];
  int *indices[2*dim];
};

#endif