amroc/amr/AMRFixup.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus

#ifndef AMROC_AMRFIXUP_H
#define AMROC_AMRFIXUP_H

#include "Fixup1.h"
#include "Fixup2.h"
#include "Fixup3.h"

#include "Timing.h"

#include <string>

template <class VectorType, class FixupType, int dim>
class AMRFixup : public Fixup<VectorType,FixupType,dim>,
    public FixupOps<VectorType,FixupType,dim> {

  typedef Fixup<VectorType,FixupType,dim> base;
  typedef FixupOps<VectorType,FixupType,dim> ops;
  typedef typename VectorType::InternalDataType DataType;

public:
  typedef GridData<VectorType,dim> vec_grid_data_type;

  typedef GridData<VectorType,minus_1<dim>::dim> ld_vec_grid_data_type;
  typedef GridData<FixupType,minus_1<dim>::dim>  ld_fixup_grid_data_type;

  typedef GridFunction<VectorType,dim> vec_grid_fct_type;  
  typedef GridFunction<FixupType,minus_1<dim>::dim> ld_fixup_grid_fct_type;  

  AMRFixup() : base() {
    _u = (vec_grid_fct_type *) 0;
    sprintf(GFName,"u-fixup");
  }

  virtual ~AMRFixup() {}

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"AMRFixup");
    RegisterAt(base::LocCtrl,"Equations",base::_FixupEquations);
  } 
  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }

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

    if (!U().has_children(Time,Level,c)) return;
    DCoords dx = base::GH().worldStep(U().GF_StepSize(Level));
    const int refinefactor = RefineFactor(base::GH(),Level);

    int d;
    DataType Fac = -dt;
    for (d=0; d<dim; d++) 
      Fac *= dx(d);
    
#ifdef DEBUG_PRINT_WHOLE_FIXUP
    ( comm_service::log() << "\n*** Saving Coarse Fluxes - " 
      << U().interiorbbox(Time,Level,c) << " ***\n").flush();
#endif
    for (register int j=0; j<U().children(Time,Level,c); j++) {
      BBox work(U().childbox(Time,Level,c,j)*U().interiorbbox(Time,Level,c));
      if (!work.empty()) {
        const int& idx = U().childidx(Time,Level,c,j);
        for (d=0; d<2*dim; d++) 
          if ((mdim==0 || d/2==mdim-1) && ValidSide(Time,Level+1,idx,d)) {
            BBox bb(FixupBBox(coarsen(U().interiorbbox(Time,Level+1,idx),refinefactor),d)*work);
            if (!bb.empty()) {
              bb = FluxBBox(bb,d);
#ifdef DEBUG_PRINT_WHOLE_FIXUP
              ( comm_service::log() << "To " << U().interiorbbox(Time,Level+1,idx) 
                << " Side: " << d << " using " << bb << " into " 
                << base::F(d)(Time,Level+1,idx).bbox() << "\n").flush();
#endif
              (*flux[d]).multiply(Fac/dx(d/2), bb);
              copyf_to(base::F(d)(Time,Level+1,idx),(*flux[d]),bb,d);
              (*flux[d]).divide(Fac/dx(d/2), bb);
#ifdef DEBUG_PRINT_WHOLE_FIXUP
              ops::AlignBBox(bb,d);
              base::debug_print_ld(base::F(d)(Time,Level+1,idx),bb); 
              ( comm_service::log() << "\n" ).flush();
#endif
            }
          }
      }
    }
  }

  virtual void AddFluxes(const int Time, const int Level, const int c,
                         vec_grid_data_type* flux[], 
                         const double tc, const double tf, 
                         const double dtc, const double dtf, 
                         const int& mdim) {
    AddIntercellFluxes(Time, Level, c, flux, dtf, mdim);
  }
  
  virtual void AddIntercellFluxes(const int Time, const int Level, const int c,
                                  vec_grid_data_type* flux[], const double dt,
                                  const int& mdim) {
    
    DCoords dx = base::GH().worldStep(U().GF_StepSize(Level));

    int d;
    DataType Fac = dt;
    for (d=0; d<dim; d++) 
      Fac *= dx(d);

#ifdef DEBUG_PRINT_WHOLE_FIXUP
    ( comm_service::log() << "\n*** Adding Higher Order Fluxes - " 
      << U().interiorbbox(Time,Level,c) << " ***\n").flush();
#endif
    for (d=0; d<2*dim; d++) 
      if ((mdim==0 || d/2==mdim-1) && ValidSide(Time,Level,c,d)) {
        BBox bb(FluxBBox(U().interiorbbox(Time,Level,c), d));
        (*flux[d]).multiply(Fac/dx(d/2), bb);
        addf_to(base::F(d)(Time,Level,c),(*flux[d]),bb,d);
        (*flux[d]).divide(Fac/dx(d/2), bb);
#ifdef DEBUG_PRINT_WHOLE_FIXUP
        base::debug_print_ld(base::F(d)(Time,Level,c));
        ( comm_service::log() << "\n").flush();
#endif
      }
  }
  
  virtual void Correction(const int Time, const int WTime, const int Level, const double t, const double dt) {
    int t_sten = 0;   
    int s_sten = 1;       
    vec_grid_fct_type uh(GFName, t_sten, s_sten, CurrentTime(base::GH(),Level), Level, 
                         base::GH(), DAGHCellCentered, 0, 1, DAGHCommSimple, 
                         DAGHNoBoundary, DAGHNoAdaptBoundary, DAGHNoExternalGhost);

#ifdef DEBUG_PRINT_WHOLE_FIXUP
    ( comm_service::log() << "\n*********** Conservative Fixup ***********\n").flush();
#endif

    for (int d=0; d<2*dim; d++) 
      ParallelFixup(Time, Level, d, uh);
  }  

  inline void SetGridFunctions(vec_grid_fct_type* u) { _u = u; }
  inline vec_grid_fct_type& U() { return *_u; }
  inline const vec_grid_fct_type& U() const { return *_u; }

protected:  
  int ValidSide(const int Time, const int Level, const int c, const int s) const {
    return (U().has_adaptiveboundary(Time,Level,c,s));
  }
    
  void ParallelFixup(const int Time, const int Level, 
                     const int d, vec_grid_fct_type& uh) {

    int TimeF = CurrentTime(base::GH(),Level+1);
    int TimeC = CurrentTime(base::GH(),Level);
    const int refinefactor = RefineFactor(base::GH(),Level);
    DCoords dx = base::GH().worldStep(U().GF_StepSize(Level));

    DataType Fac = 1.0;
    for (int dd=0; dd<dim; dd++) 
      Fac *= dx(dd);

    forall (uh,TimeC,Level,c)  
      uh(TimeC,Level,c) = 0.0;
    end_forall

    forall (U(),TimeF,Level+1,c) 
      if (ValidSide(TimeF,Level+1,c,d)) {
        BBox to(FixupBBox(coarsen(U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
        for (register int j=0; j<U().parents(TimeF,Level+1,c) ; j++) {
          BBox wto(coarsen(U().parentbox(TimeF,Level+1,c,j),refinefactor)*to);
          const int& idx = U().parentidx(TimeF,Level+1,c,j);
          wto *= uh(TimeC,Level,idx).bbox();
          if (!wto.empty()) {
            copyf_from(uh(TimeC,Level,idx),wto,base::F(d)(TimeF,Level+1,c),d);
            uh(TimeC,Level,idx).multiply((d%2==0 ? -1.0 : 1.0) / Fac,wto);
#ifdef DEBUG_PRINT_WHOLE_FIXUP
            ( comm_service::log() << "Writing at Side: " << d << " of " 
              << U().interiorbbox(TimeF,Level+1,c) << " using " << wto << "\n").flush();
            base::debug_print(uh(TimeC,Level,idx), wto);
            ( comm_service::log() << "\n" ).flush(); 
#endif  
          }
        }
      }
    end_forall 

    END_WATCH_FIXUP_WHOLE
    START_WATCH
      Sync(uh,TimeC,Level);
    END_WATCH_FIXUP_SYNC
    START_WATCH

    forall (U(),Time,Level,c)  
      BBox to = U().interiorbbox(Time,Level,c);
      BBox from(to);
      from.shift(d/2, d%2 ? -1 : 1);
      U()(Time,Level,c).plus(uh(TimeC,Level,c), to, from);
    end_forall 
  }

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

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

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

protected:
  inline void copyf_to(ld_fixup_grid_data_type &target, const BBox &to, 
                       const vec_grid_data_type &source, 
                       const BBox &from, const int s) 
    { ops::copy_to(target, to, source, from, s); }
  inline void copyf_to(ld_fixup_grid_data_type &target, 
                       const vec_grid_data_type &source, 
                       const BBox &where, const int s) 
    { ops::copy_to(target, source, where, s); }

  inline void copyf_from(vec_grid_data_type &target, const BBox &to, 
                         const ld_fixup_grid_data_type &source, 
                         const BBox &from, const int s) 
    { ops::copy_from(target, to, source, from, s); }
  inline void copyf_from(vec_grid_data_type &target, const BBox &where, 
                         const ld_fixup_grid_data_type &source, const int s) 
    { ops::copy_from(target, where, source, s); }

  inline void addf_to(ld_fixup_grid_data_type &target, const BBox &to, 
                      const vec_grid_data_type &source, 
                      const BBox &from, const int s) 
    { ops::add_to(target, to, source, from, s); }
  inline void addf_to(ld_fixup_grid_data_type &target, 
                      const vec_grid_data_type &source, 
                      const BBox &where, const int s) 
    { ops::add_to(target, source, where, s); }

  inline void addf_from(vec_grid_data_type &target, const BBox &to, 
                         const ld_fixup_grid_data_type &source, 
                         const BBox &from, const int s) 
    { ops::add_from(target, to, source, from, s); }
  inline void addf_from(vec_grid_data_type &target, const BBox &where, 
                         const ld_fixup_grid_data_type &source, const int s) 
    { ops::add_from(target, where, source, s); }

  inline void addf_to(ld_fixup_grid_data_type &target, const BBox &to, 
                       const ld_vec_grid_data_type &source, 
                       const BBox &from) 
    { ops::add_to(target, to, source, from); }  
  inline void addf_to(ld_fixup_grid_data_type &target, 
                       const ld_vec_grid_data_type &source, 
                       const BBox &where) 
    { ops::add_to(target, source, where); }  

protected:
  vec_grid_fct_type* _u;
  char GFName[DAGHBktGFNameWidth];
};

#endif

---