AMRSolver.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_AMRSOLVER_H
#define AMROC_AMRSOLVER_H

#include "AMRSolverBase.h"
#include "DAGHCluster.h"

#include <iostream>
#include <sstream>
#include <string>
#include <cstdio>
#include <cstdlib>

#define AMRSolverRegridFalse             (0)
#define VERY_LARGE_CFL                 (1.5)
#define MAXOVERLAPWIDTH                (100) 

//******************************************************************************
// flag_cells_by_error(), flag_cells_by_difference() and the fixup make use of
// u(Time+TStep). Ensure that the next time step is only used when all data,
// especially the adaptive boundary conditions (-> RecomposeHierarchy() ) have 
// been set.
// If problems occur, the use of a working function might be unavoidable.
//******************************************************************************

template <class VectorType, class FixupType, class FlagType, int dim>
class AMRSolver : public AMRSolverBase<VectorType,FixupType,FlagType,dim> {
  typedef AMRSolverBase<VectorType,FixupType,FlagType,dim> base;

  typedef Vector<FlagType,1> VFlagType;
  typedef GridData<VFlagType,dim> vflag_grid_data_type;

public:
  typedef GridData<FlagType,dim>  flag_data_type;
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::vec_grid_data_type vec_grid_data_type;  
  typedef typename base::grid_fct_type grid_fct_type;
  typedef typename base::grid_data_type grid_data_type;

  typedef typename base::integrator_type integrator_type;
  typedef typename base::initial_condition_type initial_condition_type;
  typedef typename base::boundary_conditions_type boundary_conditions_type;
  typedef typename base::leveltransfer_type leveltransfer_type;
  typedef typename base::flagging_type flagging_type;
  typedef typename base::fixup_type fixup_type;

  AMRSolver(integrator_type& integ, initial_condition_type& init,
            boundary_conditions_type& bc) : base(integ, init, bc), 
      FixupPar(1), RegridEvery(1), MinEfficiency(0.7), 
      BlockWidth(4), OverlapWidth(0), BufferWidth(2), NestingBuffer(1), NoTimeRefine(0),
      PlotFlags(0) {

    t = (double *)NULL; dt = (double *)NULL; cfl_new = (double *)NULL; 
    AllocError::SetTexts("AMRSolver","Constructor");
  }

  virtual ~AMRSolver() {    
    if (t) delete [] t;
    if (dt) delete [] dt;
    if (cfl_new) delete [] cfl_new;
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"AMRSolver");
    RegisterAt(base::LocCtrl,"DoFixup",FixupPar);
    RegisterAt(base::LocCtrl,"NoTimeRefine",NoTimeRefine);
    RegisterAt(base::LocCtrl,"RegridEvery",RegridEvery);
    RegisterAt(base::LocCtrl,"MinEfficiency",MinEfficiency);
    RegisterAt(base::LocCtrl,"BufferWidth",BufferWidth);
    RegisterAt(base::LocCtrl,"BlockWidth",BlockWidth);
    RegisterAt(base::LocCtrl,"OverlapWidth",OverlapWidth);    
    RegisterAt(base::LocCtrl,"NestingBuffer",NestingBuffer);
    RegisterAt(base::LocCtrl,"OutputFlags",PlotFlags);
    base::register_at(base::LocCtrl,"");
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual void SetupData() {
    base::SetupData();

    if (base::_Flagging && MaxLevel(base::GH())>=1) 
      base::Flagging_().SetBufferwidth(std::min(std::max(BufferWidth,OverlapWidth-1),
                                                MAXOVERLAPWIDTH));
    else
      RegridEvery = 0; 
    if (FixupPar == 0) 
      base::SetFixup((fixup_type*) 0);
      
    t = new double[MaxLevel(base::GH())+1];
    dt = new double[MaxLevel(base::GH())+1];
    cfl_new = new double[MaxLevel(base::GH())+1];

    if (FixupPar && NestingBuffer<1 && RegridEvery!=AMRSolverRegridFalse) {
      if (MY_PROC == VizServer) 
        std::cout << "   Fixup requires a NestingBuffer>0 ! Using Default." << std::endl;
      NestingBuffer = 1;
    }
    base::GH().DAGH_SetNestingBuffer(NestingBuffer);

    if (base::_Flagging) 
      base::_Flagging->SetGridFunctions(base::_u, base::_u_sh, base::_work, base::_work_sh);
    if (base::_Fixup) 
      base::_Fixup->SetGridFunctions(base::_u);
  }

  virtual void Initialize_(const double& dt_start) {   
    START_WATCH

    t[0] = 0.0; dt[0] = dt_start; cfl_new[0] = 0.0;
    for (register int l=1; l<=MaxLevel(base::GH()); l++) {
      t[l] = 0.0;
      cfl_new[l] = 0.0;
    }

    base::U().GF_DeleteStorage(1);
    if (base::ErrorEstimation()) 
      base::Ush().GF_DeleteStorage(1);

    //---------------------------------------------------------------
    // Calculate first grid distribution  
    //---------------------------------------------------------------
    int lev = 0;
    if (RegridEvery != AMRSolverRegridFalse) {
      while (lev<=FineLevel(base::GH()) && lev<MaxLevel(base::GH())) {
        //------------------------------------------------------------------
        // Set up initial data  
        //------------------------------------------------------------------                      
#ifdef DEBUG_PRINT_AMRSOLVER
        ( comm_service::log() <<  "  *** Setting flags at Level: " << lev << "\n" ).flush();
#endif
        //------------------------------------------------------------------
        // Initializing all levels ensures correct prolongation
        //------------------------------------------------------------------                      
        for (register int l=0; l<=FineLevel(base::GH()); l++) {
          base::InitialCondition_().Set(base::U(), base::Work(), l);
          SetBndry(base::U(),0,lev,0.0);
          if (base::ErrorEstimation()) {
            base::InitialCondition_().Set(base::Ush(), base::Worksh(), l);
            SetBndry(base::Ush(),0,lev,0.0);
          }
        }
        BBoxList bblist;
        RegridLevel(0, 0, bblist, false);
 
        END_WATCH_INITIALIZATION
        START_WATCH
          RecomposeHierarchy(base::GH());
        END_WATCH_RECOMPOSING_WHOLE
        START_WATCH       
        lev++;
      }

      if (FineLevel(base::GH())>0) {
        for (register int l=0; l<=FineLevel(base::GH()); l++) {
          base::InitialCondition_().Set(base::U(), base::Work(), l);
          SetBndry(base::U(),0,l,0.0);
          if (base::ErrorEstimation()) {
            base::InitialCondition_().Set(base::Ush(), base::Worksh(), l);
            SetBndry(base::Ush(),0,l,0.0);
          }
        }
        BBoxList bblist;
        forall (base::U(),0,FineLevel(base::GH()),c)      
          bblist.add(coarsen(base::U().interiorbbox(0,FineLevel(base::GH()),c),
                             RefineFactor(base::GH(),FineLevel(base::GH())-1)));
        end_forall
        RegridLevel(0, 0, bblist, true);
        END_WATCH_INITIALIZATION
        START_WATCH
          RecomposeHierarchy(base::GH());
        END_WATCH_RECOMPOSING_WHOLE
        START_WATCH
      }
    }
 
    base::U().GF_CreateStorage(1);
    if (base::ErrorEstimation()) 
      base::Ush().GF_CreateStorage(1);
    
    for (lev=0; lev<=FineLevel(base::GH()); lev++) { 
#ifdef DEBUG_PRINT_AMRSOLVER
      ( comm_service::log() <<  "  *** Initializing Level: " << lev << " *** \n" ).flush();
#endif 
      t[lev] = 0.0;
      base::InitialCondition_().Set(base::U(), base::Work(), lev);
      SetBndry(base::U(),0,lev,t[lev]);
      if (base::ErrorEstimation()) {
        base::InitialCondition_().Set(base::Ush(), base::Worksh(), lev);
        SetBndry(base::Ush(),0,lev,t[lev]);     
      }
    }    

    SetTimeSpecs(base::GH(), 0.0, dt_start, 2);

    for (lev=FineLevel(base::GH()); lev>=1; lev--) {    
      Restrict(base::U(), 0, lev, 0, lev-1);
      if (base::ErrorEstimation()) 
        Restrict(base::Ush(), 0, lev, 0, lev-1);
    }

    END_WATCH_INITIALIZATION
  }

  virtual bool Restart_(const char* CheckpointFile) {
    START_WATCH
      bool ret = RecomposeHierarchy(base::GH(), CheckpointFile);
      if (ret) {
        double NextTime;
        for (register int l=0; l<=FineLevel(base::GH()); l++) {
          base::GH().DAGH_GetTimeSpecs(t[l], NextTime);
          int Time = CurrentTime(base::GH(),l);
          SetPhysicalTime(base::U(),Time,l,t[l]);
          if (base::ErrorEstimation()) SetPhysicalTime(base::Ush(),Time,l,t[l]);
        }
      }
    END_WATCH_INITIALIZATION
      return ret;
  }

  virtual void Checkpointing_(const char* CheckpointFile) {
    Checkpoint(base::GH(), CheckpointFile);   
  }

  virtual double Tick(int VariableTimeStepping, const double dtv[], 
                      const double cflv[], int& Rejections) {

    double told = t[0];
    double cfl_max;
    Rejections = 0;
    do {
      std::stringstream* RestartStorage = (std::stringstream *)NULL;
      if (VariableTimeStepping==AutomaticTimeStepsRestart) {
        AllocError::SetTexts("AMRSolver","Memory for restart");      
        RestartStorage = new std::stringstream;
        base::GH().DAGH_Checkpoint(*RestartStorage); 
      }  

#ifdef DEBUG_PRINT_AMRSOLVER
      double Current_Time, Next_Time;
      base::GH().DAGH_GetTimeSpecs(Current_Time, Next_Time);
      ( comm_service::log() <<  " --- Iteration: " << StepsTaken(base::GH(),0)+1 << " ---"
        << "   t = " << Current_Time << "   dt = " 
        << DeltaT(base::GH(), 0) << "\n" ).flush();
#endif
#ifdef DEBUG_PROFILE
      ( std::cout << "\n" ).flush();
#endif
      
      for (register int l=1; l<=MaxLevel(base::GH()); l++)
        t[l] = t[0];

      AdvanceLevel(0, RegridEvery, 
                   !(RegridEvery != AMRSolverRegridFalse ? 
                     (StepsTaken(base::GH(),0)%RegridEvery==0 && StepsTaken(base::GH(),0)>0) : true),
                   base::ErrorEstimation() && (RegridEvery != AMRSolverRegridFalse ? 
                     (StepsTaken(base::GH(),0)+1)%RegridEvery==0 : false), FixupPar > 0,
                   !(base::RedistributeEvery>0 && StepsTaken(base::GH(),0)%base::RedistributeEvery!=0),
                   base::RedistributeEvery<=0); 

      cfl_max = 0.0;
      for (register int l=0; l<=FineLevel(base::GH()); l++)
        cfl_max = cfl_max > cfl_new[l] ? cfl_max : cfl_new[l];
      
#ifdef DAGH_NO_MPI
#else
      if (comm_service::dce() && comm_service::proc_num() > 1) { 
        
        int num = comm_service::proc_num();
        
        int R;
        int size = sizeof(double);
        void *values = (void *) new char[size*num];
        
        R = MPI_Allgather(&cfl_max, size, MPI_BYTE, values, size, MPI_BYTE, 
                          comm_service::comm(base::U().GF_Id()));
        if ( MPI_SUCCESS != R ) 
          comm_service::error_die( "AMRSolver::cfl_max", "MPI_Allgather", R );

        for (int i=0;i<num;i++) {
          double tmax = *((double *)values + i);
          cfl_max = cfl_max > tmax ? cfl_max : tmax;
        }
        if (values) delete [] (char *)values;
      }
#endif
     
#ifdef DEBUG_PRINT_AMRSOLVER
      ( comm_service::log() << "   Levels = " << FineLevel(base::GH())+1 
        << "   max. CFL = " << cfl_max << "\n" ).flush();
#endif
#ifdef DEBUG_PROFILE
      ( std::cout << "\n" ).flush();
#endif

      if (cfl_max > cflv[0]) 
        if (VariableTimeStepping==AutomaticTimeStepsRestart) {
          if (RestartStorage!=(std::stringstream *)NULL) {          
            base::GH().DAGH_RecomposeHierarchy(*RestartStorage);
            for (int lev=0; lev<=FineLevel(base::GH()); lev++)  
              t[lev] = told;
            Rejections++;
          }
          else
            std::cerr << "   Restart failure!" << std::endl << std::flush ;         
          if (cfl_max > VERY_LARGE_CFL)
            std::cerr << "   CFL>" << (double)VERY_LARGE_CFL << "  " << std::flush;
        }
        else if (cfl_max > 1.0)
          std::cerr << "   CFL>1.0  " << std::flush;
        
      if (VariableTimeStepping!=FixedTimeSteps) 
        if (cfl_max == 0.0)
          SetTimeSpecs(base::GH(), t[0], t[0]+dtv[1], 2);
        else {
          double dt_predict = (DeltaT(base::GH(), 0) / cfl_max) * cflv[1];
          SetTimeSpecs(base::GH(), t[0], t[0] + 
                       (dtv[1] < dt_predict ? dtv[1] : dt_predict), 2);
        }
      else
        SetTimeSpecs(base::GH(), t[0], t[0]+dtv[0], 2);

      if (RestartStorage != (std::stringstream *)NULL)
        delete RestartStorage;

    } while (t[0] <= told);

    return cfl_max;
  }  

  virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level, 
                                double t, double dt, bool DoFixup, double tc, 
                                const int which) {

    double cfl = 0, cfl_step;
    int mpass = 0; 
    int TStep = TimeStep(u,Level);
    if (base::Integrator_().MaxIntegratorPasses() > 0) mpass = 1;
    do {

      cfl_step = UpdateLevel(u,Time,Level,t,dt,DoFixup,tc,which,mpass);
      cfl = (cfl > cfl_step ? cfl : cfl_step);

      // Set boundary values again at t=Time. As the partially updated time
      // step is stored at Time+TStep the two time steps have to swaped. 
      if (mpass>0 && mpass<base::Integrator_().MaxIntegratorPasses()) {
        SwapTimeLevels(u,Level,Time,Time+TStep);
        SetBndry(u,Time,Level,t);
        SwapTimeLevels(u,Level,Time,Time+TStep);
      }

      mpass++;
    } while (mpass <= base::Integrator_().MaxIntegratorPasses());

#ifdef DEBUG_PRINT_AMRSOLVER
    ( comm_service::log() << "IntegrateLevel(): Level" << Level << "   local CFL = " << cfl << "\n" ).flush();
#endif
    return cfl;
  }

  //---------------------------------------------------------------
  // Take one numerical update step on u. u(Time, Level) remains unchanged. 
  //---------------------------------------------------------------

  virtual double UpdateLevel(vec_grid_fct_type& u, const int& Time, const int& Level, 
                             double t, double dt, bool DoFixup, double tc, 
                             const int& which, const int& mpass) {

    char errtext[256], whichtext[16];
    if (which<=0 || which>2) std::sprintf(whichtext,"Main");
    else if (which==1) std::sprintf(whichtext,"Estimate");
    else std::sprintf(whichtext,"Shadow");
    std::sprintf(errtext,"Flux allocation for %s",whichtext);
    AllocError::SetTexts("AMRSolver<dim>",errtext);      
    base::Integrator_().SetGridFunction(&u);

    int TStep = TimeStep(u,Level);
    double cfl = 0, cfl_grid;

    forall (u,Time,Level,c)      
      vec_grid_data_type& u_old = u(Time, Level, c);
      vec_grid_data_type& u_new = u(Time+TStep, Level, c); 

      vec_grid_data_type** f;
      base::Integrator_().AllocGridFluxes(u_old.bbox(), f);

#ifdef DEBUG_PRINT_INTEGRATOR
      if (mpass<=1) ( comm_service::log() << "      " << which << ": " ).flush();
#endif
      START_WATCH
        cfl_grid = base::Integrator_().CalculateGrid(u_new, u_old, f, Level, t, dt, mpass);
      END_WATCH_INTEGRATION(which)  
#ifdef DEBUG_PROFILE
      if (mpass<=1) ( std::cout << "." ).flush();
#endif
 
      cfl = (cfl > cfl_grid ? cfl : cfl_grid);

      if (Level>0 && DoFixup && cfl_grid>0 && cfl_grid<=VERY_LARGE_CFL) {
        // add_first_order_fluxes() solves a Riemann-problem on the basis of
        // u(CurrentTime(GH,Level-1),Level-1) at fine-coarse interfaces. 
        // Do not change u(CurrentTime(GH,Level-1),Level-1)!
        START_WATCH
          base::Fixup_().AddFluxes(Time, Level, c, f, tc, t, dt, mpass);
        END_WATCH(FIXUP_WHOLE)
      }
          
      if (Level<FineLevel(base::GH()) && DoFixup) {
        START_WATCH
          // The negation here is necessary to capture NaN's!
          if (!(cfl_grid>0 && cfl_grid<=VERY_LARGE_CFL)) 
            base::Integrator_().ResetGridFluxes(f);
          base::Fixup_().SaveFluxes(Time, Level, c, f , dt, mpass);
        END_WATCH(FIXUP_WHOLE)
      }

      base::Integrator_().DeAllocGridFluxes(f);
    end_forall

#ifdef DEBUG_PRINT_AMRSOLVER
    ( comm_service::log() << "UpdateLevel(): Level" << Level 
      << "  mpass = " << mpass << "   local CFL = " << cfl << "\n" ).flush();
#endif
    return cfl;
  }

  //---------------------------------------------------------------
  // Set boundary conditions
  //
  // Interpolation at fine-coarse interfaces is done in the following way:
  // k=0: Space prolongation from 
  //      u(CurrentTime(GH,Level-1),Level-1)
  // k>0: Space prolongation of 
  //      u(CurrentTime(GH,Level-1)+TimeStep(GH,Level-1),
  //        Level-1) and (CurrentTime(GH,Level-1,) and
  //      time interpolation
  // Ghostcells are filled in the following order:
  // 1. forall (u,Time,Level,c) 
  //      Interpolate completely into all sides abutting a patch of the coarser level.
  // 2. forall (u,Time,Level,c) 
  //      Override previous values in ghostcells overlying a patch of the same level.
  // 3. forall (u,Time,Level,c)
  //      Set physical boundary conditions at sides outside of the domain.
  // BoundaryUpdate() calls all necessary functions in the right order!
  // RecomposeHierarchy() calls BoundaryUpdate() on all available time steps
  // and levels. Therefore, we can skip the setting of boundary conditions on fine 
  // grid levels, if RecomposeHierarchy() has been called on the coarse base level.
  //---------------------------------------------------------------
  
  virtual void SetBndry(vec_grid_fct_type& u, const int Time, 
                        const int Level, double t) {    
    START_WATCH
      AdaptiveBoundaryUpdate(u,Time,Level);
    END_WATCH(BOUNDARIES_INTERPOLATION) 
    START_WATCH
      Sync(u,Time,Level);
    END_WATCH(BOUNDARIES_SYNC)
    START_WATCH
      ExternalBoundaryUpdate(u,Time,Level);
    END_WATCH(BOUNDARIES_EXTERNAL)
#ifdef DEBUG_AMR
    CheckLevel(u,Time,Level,t,"AMRSolver<dim>.SetBndry::after SetBndry");
#endif
  }

public:
  inline const int& Bufferwidth() const { return BufferWidth; }
  inline int Bufferwidth() { return BufferWidth; }

protected:
  virtual void RegridLevel(const int Time, const int BaseLevel, 
                           BBoxList& bblist, bool TakeListOnFinestLevel) {

#ifdef DEBUG_PRINT_AMRSOLVER
    ( comm_service::log() <<  "  *** Regridding at Level: " << BaseLevel 
      << "   t: " << Time << " Max: " << MaxLevel(base::GH()) << " fine: " 
      << FineLevel(base::GH()) << "\n" ).flush();
#endif

    int l;
    int flev = FineLevel(base::GH());
    if (flev == 0) l = 0;
    else if (flev >= MaxLevel(base::GH())) l = MaxLevel(base::GH())-1;      
    else l = flev;
    if (TakeListOnFinestLevel) {
      Refine(base::GH(),bblist,l);      
      l--;
    } 

    for (;l>=BaseLevel;l--) { 
      //---------------------------------------------------------------
      // Flag cells for refinement by error estimation
      //---------------------------------------------------------------
      START_WATCH

#ifdef DEBUG_AMR
        CheckLevel(base::U(),Time,l,t[l],"AMRSolver.RegridLevel::before Flagging");
#endif
        base::Flagging_().SetFlags(Time, l, t[l], dt[l]);

        // Ensure nesting of higher level patches. Then syncronize the flags.
        forall (base::Flagging_().Flags(),Time,l,c)  
          for (BBox *_b = bblist.first();_b;_b=bblist.next()) 
            if (!_b->empty()) {
              BBox wb = base::Flagging_().Flags().interiorbbox(Time,l,c);
              if(wb.stepsize(0)>_b->stepsize(0))
                _b->coarsen(wb.stepsize(0)/_b->stepsize(0));
              base::Flagging_().Flags()(Time,l,c).plus(200,*_b);
            }
        end_forall

        // The flags have to be synced to create the bufferzone around flagged cells
        // correctly. Without syncing no bufferzone into neighbouring patches on other
        // processors would be created.
        // The radius for the syncronization stencil is set to
        // Min(Max(BufferWidth,OverlapWidth-1),MAXOVERLAPWIDTH)
        START_WATCH
          Sync(base::Flagging_().Flags(), Time, l);
        END_WATCH(BOUNDARIES_SYNC)
      END_WATCH_FLAGGING 

      if (PlotFlags && base::LastOutputTime()==Time) {
        char flagname[DAGHBktGFNameWidth+16];
        std::sprintf(flagname,"flags"); 
        forall (base::Work(),Time,l,c)  
          flag_data_type& flag = base::Flagging_().Flags()(Time,l,c);
          vflag_grid_data_type vflag_help(flag.bbox(),(VFlagType*) flag.data());
          equals_from(base::Work()(Time,l,c), vflag_help, 0);
          VFlagType* databuf;
          vflag_help.deallocate(databuf);
        end_forall  
        base::FileOutput_().WritePlain(base::Work(), flagname, Time, l, t[l]);  
      }  

      //---------------------------------------------------------------
      // Cluster at refine 
      //---------------------------------------------------------------
#ifdef DEBUG_PRINT_INTEGRATOR
      ( comm_service::log() << "  *** Clustering at Level: " << l << "\n" ).flush();
#endif
      START_WATCH        
        BBoxList bblexclude, bblcut;
        for (BBox *bbi=base::GH().wholebaselist()->first();bbi;
             bbi=base::GH().wholebaselist()->next())      
          bblexclude.add(refine(*bbi,base::GH().refinedby(l)));
        base::GH().glb_bboxlist(bblcut,l);
        bblexclude -= bblcut;

        // After syncing, the clusterer creates the bufferzone locally.
        DAGHCluster(base::Flagging_().Flags(), Time, l, 
                    BlockWidth, OverlapWidth, BufferWidth, 
                    MinEfficiency, (FlagType)GoodPoint,
                    bblist, bblist, bblexclude);
      END_WATCH_CLUSTERING

      if (bblist.isempty()) {
        BBox empty_box;
        bblist.add(empty_box);
      }

#ifdef DEBUG_PRINT_RG_REFINE
      ( comm_service::log() << "\n************* New refine list at level " 
        << l << " *************\n" << bblist
        << "\n************* ***************** *************\n"
        ).flush();
#endif

      Refine(base::GH(),bblist,l);      
    }
  }

  virtual void RecomposeGridHierarchy(const int Time, const int Level, 
                                      bool ShadowAllowed, bool DoFixup,
                                      bool RecomposeBaseLev, bool RecomposeHighLev) {
    // Time interpolation needs two time steps only on the coarser levels.
    if (DoFixup) 
      base::Fixup_().SetMaxRecomposeLevel(Level);
    register int lev;
    for (lev=Level; lev<=MaxLevel(base::GH()); lev++) {
      base::U().GF_DeleteStorage(1,lev);
      if (ShadowAllowed) base::Ush().GF_DeleteStorage(1,lev);
    }
#ifdef DEBUG_AMR
    if (base::Integrator_().Abort()) 
      for (register int l=0; l<=FineLevel(base::GH()); l++) {
        int time = CurrentTime(base::GH(),l);
        int tstep = TimeStep(base::U(),l);
        CheckLevel(base::U(),time,l,t[l],"AMRSolver.AdvanceLevel::before Recompose");
        if (l<Level) 
          CheckLevel(base::U(),time+tstep,l,t[l],
                     "AMRSolver.AdvanceLevel::before Recompose for Time+TStep");
        if (ShadowAllowed && l<MaxLevel(base::GH())) {
          int shtstep = TimeStep(base::Ush(),l);
          CheckLevel(base::Ush(),time,l,t[l],
                     "AMRSolver.AdvanceLevel::before Recompose for Shadow");
          if (l<Level) 
            CheckLevel(base::Ush(),time+shtstep,l,t[l],
                       "AMRSolver.AdvanceLevel::before Recompose for Time+TStep for Shadow");
        }               
      }
#endif
    RecomposeHierarchy(base::GH(), (Level==0 && RecomposeBaseLev) || 
                       (Level>0 && RecomposeHighLev) ? DAGHTrue : DAGHFalse);
#ifdef DEBUG_AMR
    if (base::Integrator_().Abort()) 
      for (register int l=0; l<=FineLevel(base::GH()); l++) {
        int time = CurrentTime(base::GH(),l);
        int tstep = TimeStep(base::U(),l);
        CheckLevel(base::U(),time,l,t[l],"AMRSolver.AdvanceLevel::after Recompose");
        if (l<Level) 
          CheckLevel(base::U(),time+tstep,l,t[l],
                     "AMRSolver.AdvanceLevel::after Recompose for Time+TStep");
        if (ShadowAllowed && l<MaxLevel(base::GH())) {
          int shtstep = TimeStep(base::Ush(),l);
          CheckLevel(base::Ush(),time,l,t[l],
                     "AMRSolver.AdvanceLevel::after Recompose for Shadow");
          if (l<Level) 
            CheckLevel(base::Ush(),time+shtstep,l,t[l],
                       "AMRSolver.AdvanceLevel::after Recompose for Time+TStep for Shadow");
        }               
      }
#endif

    for (lev=Level; lev<=MaxLevel(base::GH()); lev++) {
      base::U().GF_CreateStorage(1,lev);
      if (ShadowAllowed && lev<MaxLevel(base::GH()))
        base::Ush().GF_CreateStorage(1,lev);
    }
  }

  virtual void AfterLevelStep(const int Level) {}

  virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone, 
                             bool ShadowAllowed, bool DoFixup,
                             bool RecomposeBaseLev, bool RecomposeHighLev) {
    //---------------------------------------------------------------
    // Select number of iterations to run based on current grid level
    //---------------------------------------------------------------
    int NoIterations = 0;
    if (Level==0 || NoTimeRefine==1) {
      NoIterations = 1;
      dt[Level] = DeltaT(base::GH(), 0);
    }
    else {
      NoIterations = RefineFactor(base::GH(),Level-1);
      dt[Level] = dt[Level-1] / RefineFactor(base::GH(),Level-1);
    }
    // Set dt also on higher levels to allow correct error estimation
    for (int l=Level+1; l<=MaxLevel(base::GH()); l++) 
      if (NoTimeRefine==1) dt[l]=dt[l-1];
      else dt[l] = dt[l-1] / RefineFactor(base::GH(),l-1);

    int Time;
    int TStep = TimeStep(base::U(),Level);
    double cfl;
    cfl_new[Level] = 0.0;

    for (int k=0; k<NoIterations; k++) {      
      Time = CurrentTime(base::GH(),Level);

#ifdef DEBUG_PRINT_AMRSOLVER
      ( comm_service::log() << "AdvanceLevel(): Level" << Level 
        << "  Time = " << Time << "   Iteration = " << k << "\n" ).flush();
#endif

      SetPhysicalTime(base::U(),Time+TStep,Level,t[Level]+dt[Level]);
      if (ShadowAllowed) 
        SetPhysicalTime(base::Ush(),Time+TimeStep(base::Ush(),Level),Level, 
                        t[Level]+dt[Level]*base::Ush().factor());

      bool DoRegrid = (RegridEvery!= AMRSolverRegridFalse ?
                       Level<MaxLevel(base::GH()) && k%RegridEvery==0 && 
                       !RegridDone && t[Level]>0.0 : false);    
     
      if (Level==0 || !RegridDone) 
        SetBndry(base::U(),Time,Level,t[Level]);

      //---------------------------------------------------------------
      // regridding...
      //---------------------------------------------------------------
      if (DoRegrid) {   
        //---------------------------------------------------------------
        // Regrid all levels > Level, Level itself stays fixed
        //---------------------------------------------------------------       
        BBoxList bblist;
        RegridLevel(Time, Level, bblist, false);
        
        //---------------------------------------------------------------
        // Reorganization of hierarchy
        //---------------------------------------------------------------       
        START_WATCH
          RecomposeGridHierarchy(Time,Level,ShadowAllowed,DoFixup, 
                                 RecomposeBaseLev,RecomposeHighLev);
        END_WATCH_RECOMPOSING_WHOLE
        RegridDone = true;
      } 
        
      //---------------------------------------------------------------
      // Take one step on U(). U()(Time, Level) remains unchanged. 
      // In principle, it would be possible to skip this step after error estimation on the 
      // unchanged level. But this is not the case here, because numerical fluxes are only 
      // available for single grids.  As fine grids might change due to error estimation,
      // the saving of coarse grid fluxes would require the storage of numerical fluxes
      // on the entire coarse level.
      //---------------------------------------------------------------      
      cfl = IntegrateLevel(base::U(), Time, Level, t[Level], dt[Level], 
                           DoFixup, (Level>0 ? t[Level-1] : 0.0), 0); 
      cfl_new[Level] = cfl_new[Level] > cfl ? cfl_new[Level] : cfl;

      //---------------------------------------------------------------
      // Take one step on Ush().
      // If restriction into Ush() is done AFTER the step on
      // U(), U()(Time,Level) has to remain unchanged in a fractional 
      // step method!
      //---------------------------------------------------------------
                
      if (RegridEvery != AMRSolverRegridFalse ?
          Level<MaxLevel(base::GH()) && ShadowAllowed &&
          (Level>0 ? (k+1)%RegridEvery==0 : true) : false) {

        double dt_sh = dt[Level] * base::Ush().factor();
        SetPhysicalTime(base::Ush(),Time,Level,t[Level]);
        
        //------------------------------------------------------------------------
        // Restrict from U() into Ush()
        //------------------------------------------------------------------------
        forall (base::U(),Time,Level,cm)   
          vec_grid_data_type& u_main = base::U()(Time, Level, cm);
          forall (base::Ush(),Time,Level,cs)    
            vec_grid_data_type& u_shadow = base::Ush()(Time, Level, cs);
            BBox bb = u_shadow.bbox()*u_main.bbox();
            if (!bb.empty())
              base::LevelTransfer_().Restrict(u_main,Level,u_shadow,Level,bb);
          end_forall
        end_forall

        SetBndry(base::Ush(),Time,Level,t[Level]);
        IntegrateLevel(base::Ush(), Time, Level, t[Level], dt_sh, false, 0.0, 2); 
        
      } // Ends if statement for Ush()
        
      //---------------------------------------------------------------
      // If we are not at the finest level then go to next level
      //---------------------------------------------------------------
      
#ifdef DEBUG_PRINT_FIXUP
      VectorType mass_old;
      int d;
      DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
      if (Level == 0) {
        mass_old = Sum(base::U(),Time+TStep,Level);
        for (d=0; d<dim; d++) mass_old *= dx(d);
      }
#endif
 
      if (Level < MaxLevel(base::GH())) {
        
        // Boundary conditions have to be applied on U()(Time+TStep,Level)
        // to ensure correct time-space interpolation for Level+1.
        SetBndry(base::U(),Time+TStep,Level,t[Level]+dt[Level]);

        // Recursive Integration of next level in hierarchy
        AdvanceLevel(Level+1, RegridEvery, RegridDone, 
                     ShadowAllowed, DoFixup, RecomposeBaseLev, 
                     RecomposeHighLev);
        
#ifdef DEBUG_AMR
        CheckLevel(base::U(),Time+TStep,Level,t[Level],
                   "AMRSolver.AdvanceLevel::before Restrict and Fixup");
#endif

        if (DoFixup) {
          START_WATCH
            base::Fixup_().Correction(Time+TStep, Time, Level);   
          END_WATCH(FIXUP_CORRECTION)
        }

        // Implementation of fixup assumes that restriction is done AFTERWARDS!
        // Only values from internal cells are used. Ghostcell values are not restricted.
        Restrict(base::U(), Time+TStep, Level+1, Time+TStep, Level);

#ifdef DEBUG_AMR
        CheckLevel(base::U(),Time+TStep,Level,t[Level],
                   "AMRSolver.AdvanceLevel::after Restrict and Fixup");
#endif
      } 
#ifdef DEBUG_PROFILE
      ( std::cout << "\n" ).flush();
#endif

#ifdef DEBUG_PRINT_FIXUP
      VectorType mass_new;
      if (Level == 0) {
        ( comm_service::log() 
          << "\n*********** Difference in state variables at level " 
          << Level << " after fixup ********\n" ).flush();
        // Sum() uses only internal cells. No synchronization on Level is therefore 
        // necessary after restriction!
        mass_new = Sum(base::U(),Time+TStep,Level);
        for (d=0; d<dim; d++) mass_new *= dx(d);
        ( comm_service::log() << mass_old << "-" 
          << mass_new << "=" << mass_old - mass_new ).flush();
        ( comm_service::log() 
          << "\n**************************************************************"
          << "***********\n" ).flush();
      }
#endif
        
      //---------------------------------------------------------------
      // Move from current time to previous.
      //---------------------------------------------------------------
      CycleTimeLevels(base::U(),Level);
      if (ShadowAllowed)
        CycleTimeLevels(base::Ush(),Level);
        
      //---------------------------------------------------------------
      //  Increment time step on current level
      //  IncrCurrentTime() moves internal time-indices.
      //  Can only be called after all time-depended operations, like
      //  time-interpolation or time-cycling.
      //---------------------------------------------------------------
      t[Level] += dt[Level]; 
      IncrCurrentTime(base::GH(),Level);
      if (NoTimeRefine==1 && Level>0) 
        for (int kk=1; kk<RefineFactor(base::GH(),Level-1); kk++) 
          IncrCurrentTime(base::GH(),Level);
      Time = CurrentTime(base::GH(),Level);
      SetPhysicalTime(base::U(),Time,Level,t[Level]);
      if (ShadowAllowed) 
          SetPhysicalTime(base::Ush(),Time,Level,t[Level]+ 
                          dt[Level]*(base::Ush().factor()-1));
      RegridDone = false;

      AfterLevelStep(Level);
    }  
  } 

  void CheckLevel(vec_grid_fct_type &u, const int& time, const int& level, 
                  const double& t, const char *text) {
    if (base::Integrator_().Abort()) {
      char wtext[1024];
      std::sprintf(wtext,"\non Level=%d at Time=%d: %s\n",level,time,text);
      forall (u,time,level,c)      
        vec_grid_data_type& uw = u(time,level,c);
        base::Integrator_().CheckGrid(uw, level, uw.bbox(), t, wtext);
      end_forall
    }
  }

  BBox BoundaryBBox(const BBox& whole, const int s) {    
    int d = s/2;
#ifdef DEBUG_PRINT
    assert (d>=0 && d<dim);
#endif
    BBox bb(whole); 
    if (s%2 == 0) 
      bb.upper(d) = bb.lower(d)+(base::NGhosts()-1)*bb.stepsize(d);
    else
      bb.lower(d) = bb.upper(d)-(base::NGhosts()-1)*bb.stepsize(d);
   return bb;
  }

  int FixupPar;
  int RegridEvery;
  double MinEfficiency;         
  int BlockWidth, OverlapWidth, BufferWidth, NestingBuffer;
  int NoTimeRefine;
  int PlotFlags;

  double *t, *dt;
  double *cfl_new;
};


#endif

---