WENOFixup.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROC_WENO_FIXUP_H
#define AMROC_WENO_FIXUP_H

#include "AMRFixup.h"

template <class VectorType, class FixupType, int dim>
class WENOFixup : public AMRFixup<VectorType,FixupType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  typedef AMRFixup<VectorType,FixupType,dim> base;

public:
  typedef typename base::vec_grid_data_type vec_grid_data_type;

  WENOFixup() : base() {}

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

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

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

#ifdef DEBUG_PRINT_WHOLE_FIXUP
    ( comm_service::log() << "\n*** Saving Coarse Fluxes - " 
      << base::U().interiorbbox(Time,Level,c) << " ***\n").flush();
#endif
    for (register int j=0; j<base::U().children(Time,Level,c); j++) {
      BBox work(base::U().childbox(Time,Level,c,j)*base::U().interiorbbox(Time,Level,c));
      if (!work.empty()) {
        const int& idx = base::U().childidx(Time,Level,c,j);
        for (d=0; d<2*dim; d++) 
          if (base::ValidSide(Time,Level+1,idx,d)) {
            BBox bb(base::FixupBBox(coarsen(base::U().interiorbbox(Time,Level+1,idx),refinefactor),d)*work);
            if (!bb.empty()) {
              bb = base::FluxBBox(bb,d);
#ifdef DEBUG_PRINT_WHOLE_FIXUP
              ( comm_service::log() << "To " << base::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);

              if (mpass<=1)
                base::copyf_to(base::F(d)(Time,Level+1,idx),(*flux[d]),bb,d);
              else
                base::addf_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
              AlignBBox(bb,d);
              base::debug_print_ld(base::F(d)(Time,Level+1,idx),bb); 
              ( comm_service::log() << "\n" ).flush();
#endif
            }
          }
      }
    }
  }

  virtual void AddIntercellFluxes(const int Time, const int Level, const int c,
                                  vec_grid_data_type* flux[], const double dt,
                                  const int& mpass) {
    
    DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
    DataType Fac = dt;
    int d;
    for (d=0; d<dim; d++) 
      Fac *= dx(d);

#ifdef DEBUG_PRINT_WHOLE_FIXUP
    ( comm_service::log() << "\n*** Adding Higher Order Fluxes - " 
      << base::U().interiorbbox(Time,Level,c) << " ***\n").flush();
#endif
    for (d=0; d<2*dim; d++) 
      if (base::ValidSide(Time,Level,c,d)) {
        BBox bb(base::FluxBBox(base::U().interiorbbox(Time,Level,c), d));
        (*flux[d]).multiply(Fac/dx(d/2), bb);
        base::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
      }
  }
  
};

#endif