ClpFixup.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROC_CLP_FIXUP_H
#define AMROC_CLP_FIXUP_H

#include "AMRFixup.h"

// -*- C++ -*-

#include <iostream>
#include <cstdlib>

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

template <class VectorType, int dim>
class ClpFixupVecOps : private FixupOps<VectorType,VectorType,dim> {
public:
  ClpFixupVecOps() {}
protected:
  inline void copyv_to(GridData<VectorType,minus_1<dim>::dim> &target, 
                       const GridData<VectorType,dim> &source, 
                       const BBox &where, const int s) 
    { FixupOps<VectorType,VectorType,dim>::copy_to(target, source, where, s); }
};

#define f_rcflx FORTRAN_NAME(rcflx, RCFLX)
extern "C" {
  void f_rcflx(const INTEGER& mflx);
}

template <class VectorType, class FixupType, class AuxVectorType, int dim>
class ClpFixup : 
  public AMRFixup<VectorType,FixupType,dim>,
  public ClpFixupVecOps<VectorType,dim> {

  typedef typename VectorType::InternalDataType DataType;
  typedef AMRFixup<VectorType,FixupType,dim> base;
  typedef ClpFixupVecOps<VectorType,dim> ops;

public:
  typedef ClpIntegrator<VectorType,AuxVectorType,dim> integrator_type; 
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::ld_vec_grid_data_type ld_vec_grid_data_type;

  ClpFixup(integrator_type& integ) : base(), _Integrator(integ) {}

protected:
  void CorrectFirstOrder(const int Time, const int Level, const int c, 
                         const double tc, const double tf,
                         const double dt, const int& mdim) {
    
    int TimeC = CurrentTime(base::GH(),Level-1);
    DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
    const int refinefactor = RefineFactor(base::GH(),Level-1);
#ifdef DEBUG_PRINT_WHOLE_FIXUP
    ( comm_service::log() << "\n*** Adding First Order Fluxes - "
      << base::U().interiorbbox(Time,Level,c) << " ***\n").flush();
#endif
      
    double Fac[2*dim]; 
    int d;
    for (d=0; d<dim; d++) {
      Fac[2*d] = dt; Fac[2*d+1] = -dt;
      for (int dd=0; dd<dim; dd++)
        if (d != dd) {
          Fac[2*d] *= dx(dd); Fac[2*d+1] *= dx(dd);
        }
    }

    for (d=0; d<2*dim; d++) 
      if ((mdim==0 || d/2==mdim-1) && base::ValidSide(Time,Level,c,d)) {
#ifdef DEBUG_PRINT_WHOLE_FIXUP
        ( comm_service::log() << "Side:  " << d << "\n").flush();
#endif    
        BBox bb(base::StateVecBBox(base::U().interiorbbox(Time,Level,c), d));
        BBox bbC(base::StateVecBBox(coarsen(base::U().interiorbbox(Time,Level,c),
                                      refinefactor),d));      
        BBox workbb(bb); base::AlignBBox(workbb, d);
        
        ld_vec_grid_data_type LeftState(workbb); 
        ld_vec_grid_data_type RightState(workbb); 
        ld_vec_grid_data_type NewState(workbb); 
        BBox bbl(bb), bbr(bb);
        double dtl = dt, dtr = dt, tl, tr;
      
        /* Collect data for Riemann problems at coarse-fine boundaries */
        if (d % 2 == 0) {
          /* Fine grid data comes from the current grid */
#ifdef DEBUG_AMR
          if (Integrator_().Abort()) 
            Integrator_().CheckGrid(base::U()(Time,Level,c), Level, bb, tf,
               "ClpFixup2::CorrectFirstOrder::LeftState,fine");
#endif
          ops::copyv_to(LeftState, base::U()(Time,Level,c), bb, d);
          tl = tf;

          /* Coarse grid data comes from all parent grids */
          for (register int j=0; j<base::U().parents(Time,Level,c) ; j++) {
            BBox wbbC(coarsen(base::U().parentbox(Time,Level,c,j),refinefactor)*bbC);
            if (!wbbC.empty()) {
#ifdef DEBUG_AMR
              if (Integrator_().Abort()) 
                Integrator_().CheckGrid(base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), Level-1, 
                                        wbbC, tc, "ClpFixup2::CorrectFirstOrder::RightState,coarse");
#endif
              ops::copyv_to(RightState, base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), wbbC, d);
            }
          }
          tr = tc;
          dtr      *= refinefactor; 
          bbr.coarsen(d/2, refinefactor);
        }
        else {
          /* Coarse grid data comes from all parent grids */
          for (register int j=0; j<base::U().parents(Time,Level,c) ; j++) {
            BBox wbbC(coarsen(base::U().parentbox(Time,Level,c,j),refinefactor)*bbC);
            if (!wbbC.empty()) {
#ifdef DEBUG_AMR
              if (Integrator_().Abort()) 
                Integrator_().CheckGrid(base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), Level-1, 
                                        wbbC, tc, "ClpFixup2::CorrectFirstOrder::LeftState,coarse");
#endif
              ops::copyv_to(LeftState, base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), wbbC, d);
            }
          }
          tl = tc;
          dtl      *= refinefactor; 
          bbl.coarsen(d/2, refinefactor);
        
          /* Fine grid data comes from the current grid */
#ifdef DEBUG_AMR
          if (Integrator_().Abort()) 
            Integrator_().CheckGrid(base::U()(Time,Level,c), Level, bb, tf,
               "ClpFixup2::CorrectFirstOrder::RightState,fine");
#endif
          ops::copyv_to(RightState, base::U()(Time,Level,c), bb, d);
          tr = tf;
        }
      
        ((integrator_type&)Integrator_()).CalculateRP(LeftState, bbl, tl, dtl,  
           RightState, bbr, tr, dtr, NewState, d/2+1);      
      
#ifdef DEBUG_PRINT_WHOLE_FIXUP
        base::debug_print_ld(base::F(d)(Time,Level,c));
        ( comm_service::log() << "\n").flush();
        ( comm_service::log() << "LeftState: \n").flush();
        base::debug_print_ldv(LeftState); ( comm_service::log() << "\n").flush();
        ( comm_service::log() << "RightState: \n").flush();
        base::debug_print_ldv(RightState); ( comm_service::log() << "\n").flush();
        ( comm_service::log() << "NewState: \n").flush();
        base::debug_print_ldv(NewState); ( comm_service::log() << "\n").flush();
#endif
        NewState.multiply(Fac[d],workbb);
        base::addf_to(base::F(d)(Time,Level,c),NewState,workbb);
#ifdef DEBUG_PRINT_WHOLE_FIXUP
        base::debug_print_ld(base::F(d)(Time,Level,c));
        ( 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 dt, const int& mdim) {
    int mflx;
    f_rcflx(mflx);
    if (mflx==0) {
      if (dim>1 && Integrator_().DimSplit()) {
        std::cerr << "\n\nThe fixup is incorrect, if dimensional splitting is used with a\n";
        std::cerr << "scheme returning flux differences. Set Method(3)>=0 instead.\n";
        std::cerr << "Aborting programm...\n" << std::flush;
        std::exit(-1);
      }
      CorrectFirstOrder(Time, Level, c, tc, tf, dt, mdim);
    }
    
    base::AddIntercellFluxes(Time, Level, c, flux, dt, mdim);
  }

  inline integrator_type& Integrator_() { return _Integrator; }
  inline const integrator_type& Integrator_() const { return _Integrator; }
  inline const int& NGhosts() const { return _Integrator.NGhosts(); }

private:
  integrator_type& _Integrator;
}; 

#endif