AMRGFMSolver.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROC_GFMSOLVER_H
#define AMROC_GFMSOLVER_H

#include "AMRSolver.h"
#include "GhostFluidMethod.h"
#include "Criteria/RefinePhi.h"
#include "Criteria/ResolvePhi.h"
#include "Criteria/UnflagPhi.h"

template <class VectorType, class FixupType, class FlagType, int dim>
class AMRGFMSolver : public AMRSolver<VectorType,FixupType,FlagType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  typedef AMRSolver<VectorType,FixupType,FlagType,dim> base;
  typedef AMRGFMSolver<VectorType,FixupType,FlagType,dim> self;
public:
  typedef GhostFluidMethod<VectorType,dim> gfm_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;

  typedef GFRecomposeSpecificFunc<self,VectorType,dim> gfm_recompose_functor_type;

  typedef GridData<bool,dim>  bool_grid_data_type;
  typedef GridFunction<bool,dim> bool_grid_fct_type;  

  AMRGFMSolver(integrator_type& integ, 
               initial_condition_type& init,
               boundary_conditions_type& bc) : base(integ, init, bc), _MaxRecomposeLevel(-1),
                                               _RecoverExterior(0) {      
    _bf = (bool_grid_fct_type *) 0;
    _bf_sh = (bool_grid_fct_type *) 0;
    _GFMRecomposeFunc = (gfm_recompose_functor_type *) 0;
    _GFM = (gfm_type**) 0;
    _nGFM = 0;
  }

  virtual ~AMRGFMSolver() {    
    if (_bf) delete _bf;
    if (_bf_sh) delete _bf_sh;
    if (_GFMRecomposeFunc) delete _GFMRecomposeFunc;
  }

  virtual void init() {
    base::Flagging_().AddCriterion(new RefinePhi<VectorType,FixupType,FlagType,dim>(*this));
    base::Flagging_().AddCriterion(new UnflagPhi<VectorType,FixupType,FlagType,dim>(*this));    
    base::Flagging_().AddCriterion(new ResolvePhi<VectorType,FixupType,FlagType,dim>(*this,false));
    base::Flagging_().AddCriterion(new ResolvePhi<VectorType,FixupType,FlagType,dim>(*this,true));
    base::init();
    for (register int n=0; n<NGFM(); n++) 
      GFM(n).init();
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
    for (register int n=0; n<NGFM(); n++) {
      char text[5];
      std::sprintf(text,"%d-",n+1);
      GFM(n).register_at(base::LocCtrl, std::string(text));
    }
    RegisterAt(base::LocCtrl,"RecoverExterior",_RecoverExterior);
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual void update() {
    base::update();
    for (register int n=0; n<NGFM(); n++) 
      GFM(n).update();
  }
  virtual void finish() {
    if (_bf) {
      delete _bf;
      _bf = (bool_grid_fct_type *) 0;
    }
    if (_bf_sh) {
      delete _bf_sh;
      _bf_sh = (bool_grid_fct_type *) 0;
    }
    base::finish();
    for (register int n=0; n<NGFM(); n++) 
      GFM(n).finish();
  }
  virtual void SetupData() {
    base::SetupGridHierarchy();
    for (register int n=0; n<NGFM(); n++) {
      GFM(n).SetErrorEstimation(base::ErrorEstimation()); 
      GFM(n).SetupData(base::PGH(), base::NGhosts());
    }

    int t_sten = 1;
    int s_sten = base::NGhosts();
    std::sprintf(BFName,"BFlags");
    _bf = new bool_grid_fct_type(BFName, t_sten, s_sten, base::GH(),
                                 DAGHCellCentered, DAGHNoComm, DAGHNoBoundary,
                                 DAGHNoAdaptBoundary, DAGHNoExternalGhost);
    SetCheckpointFlag(BF(), DAGHFalse);
    SetTimeAlias(BF(), -1, 1);
    if (base::ErrorEstimation()) {
      std::sprintf(BFNamesh,"BFlagssh");
      _bf_sh = new bool_grid_fct_type(BFNamesh, t_sten, s_sten, base::GH(),
                                      DAGHCellCentered, 0, DAGHShadowFactor, 
                                      DAGH_All, DAGHNoComm, DAGHNoBoundary,
                                      DAGHNoAdaptBoundary, DAGHNoExternalGhost);
      SetCheckpointFlag(BFsh(), DAGHFalse);
      SetTimeAlias(BFsh(), -1, 1);
    }

    base::SetupData();
    _GFMRecomposeFunc = new gfm_recompose_functor_type(this, &self::SetRecomposeBndry);
    base::U().GF_SetRecomposeFunc(_GFMRecomposeFunc);
    if (base::ErrorEstimation())  
      base::Ush().GF_SetRecomposeFunc(_GFMRecomposeFunc);
    for (register int n=0; n<NGFM(); n++) 
      GFM(n).SetGridFunctions(base::_u,base::_u_sh);
  }

  virtual void SetBndry(vec_grid_fct_type& u, const int Time, 
                        const int Level, double t) {    
    base::SetBndry(u,Time,Level,t);
    forall (BF(u),Time,Level,c)      
      BF(u)(Time,Level,c).equals(false);
    end_forall
    for (register int n=0; n<NGFM(); n++) {
      if (!GFM(n).Stationary() || Time==0) 
        GFM(n).SetLevelSet(u,Time,Level,t);     
      GFM(n).SetBndry(u,BF(u),Time,Level,t);
    }
  }

  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 = base::IntegrateLevel(u,Time,Level,t,dt,DoFixup,tc,which);

    int TStep = TimeStep(u,Level);     
    if (_RecoverExterior) {
      forall (u,Time,Level,c)           
        if (base::Dim() == 1) {
          BeginFastIndex1(u_old, u(Time,Level,c).bbox(),
                          u(Time,Level,c).data(), VectorType);
          BeginFastIndex1(u_new, u(Time+TStep,Level,c).bbox(),
                          u(Time+TStep,Level,c).data(), VectorType);
          BeginFastIndex1(bf, BF(u)(Time,Level,c).bbox(),
                          BF(u)(Time,Level,c).data(), bool);
          for_1 (n, BF(u)(Time,Level,c).bbox(), BF(u)(Time,Level,c).bbox().stepsize())
            if (FastIndex1(bf,n)) 
              FastIndex1(u_new,n) = FastIndex1(u_old,n);
          end_for
          EndFastIndex1(bf);
          EndFastIndex1(u_new);
          EndFastIndex1(u_old);
        }      
        else if (base::Dim() == 2) {
          BeginFastIndex2(u_old, u(Time,Level,c).bbox(),
                          u(Time,Level,c).data(), VectorType);
          BeginFastIndex2(u_new, u(Time+TStep,Level,c).bbox(),
                          u(Time+TStep,Level,c).data(), VectorType);
          BeginFastIndex2(bf, BF(u)(Time,Level,c).bbox(),
                          BF(u)(Time,Level,c).data(), bool);
          for_2 (n, m, BF(u)(Time,Level,c).bbox(), BF(u)(Time,Level,c).bbox().stepsize())
            if (FastIndex2(bf,n,m)) 
              FastIndex2(u_new,n,m) = FastIndex2(u_old,n,m);
          end_for
          EndFastIndex2(bf);
          EndFastIndex2(u_new);
          EndFastIndex2(u_old);
        }
        else if (base::Dim() == 3) {
          BeginFastIndex3(u_old, u(Time,Level,c).bbox(),
                          u(Time,Level,c).data(), VectorType);
          BeginFastIndex3(u_new, u(Time+TStep,Level,c).bbox(),
                          u(Time+TStep,Level,c).data(), VectorType);
          BeginFastIndex3(bf, BF(u)(Time,Level,c).bbox(),
                          BF(u)(Time,Level,c).data(), bool);
          for_3 (n, m, l, BF(u)(Time,Level,c).bbox(), BF(u)(Time,Level,c).bbox().stepsize())
            if (FastIndex3(bf,n,m,l)) 
              FastIndex3(u_new,n,m,l) = FastIndex3(u_old,n,m,l);
          end_for
          EndFastIndex3(bf);
          EndFastIndex3(u_new);
          EndFastIndex3(u_old);
        }
      end_forall
    }

    // Revert GFM boundary cells back to values from interior.
    // This allows for moving boundaries and enables the creation
    // of new refinement regions along the GFM boundary just using
    // interior cell values.
    // Swapping time levels also for BF(u) ensures correct flags
    // also in (BF(u(Time+TStep)).
    for (register int n=0; n<NGFM(); n++) 
      if (Level<MaxLevel(base::GH()))
        GFM(n).SetBndry(u,BF(u),Time,Level,t,true);

    SwapTimeLevels(u,Level,Time,Time+TStep);
    SwapTimeLevels(BF(u),Level,Time,Time+TStep);
    for (register int n=0; n<NGFM(); n++) 
      GFM(n).SetBndry(u,BF(u),Time,Level,t,true);
    SwapTimeLevels(u,Level,Time,Time+TStep);
    SwapTimeLevels(BF(u),Level,Time,Time+TStep);

    return cfl;
  }

  virtual void RecomposeGridHierarchy(const int Time, const int Level, 
                                      bool ShadowAllowed, bool DoFixup,
                                      bool RecomposeBaseLev, bool RecomposeHighLev) {

    // Recomposition of GridFunction requires LevelTransfer and a BoundaryConditions functor.
    // This is currently not available for Levelset functions.
    _MaxRecomposeLevel = DAGHNull;
    for (register int n=0; n<NGFM(); n++)
      GFM(n).SetMaxRecomposeLevel(_MaxRecomposeLevel);
    BF().GF_SetMaxRecomposeLevel(_MaxRecomposeLevel);
    if (base::ErrorEstimation())  
      BFsh().GF_SetMaxRecomposeLevel(_MaxRecomposeLevel);
    base::RecomposeGridHierarchy(Time,Level,ShadowAllowed,DoFixup,
                                 RecomposeBaseLev,RecomposeHighLev);
  }

  virtual void SetRecomposeBndry(vec_grid_fct_type& u, const int& Time, 
                                 const int& Level, const double& t) {
    if (Time>0 && CurrentTime(base::GH(),Level)==Time &&
        Level>_MaxRecomposeLevel) {
      forall (BF(u),Time,Level,c)      
        BF(u)(Time,Level,c).equals(false);
      end_forall
      for (register int n=0; n<NGFM(); n++) {
        GFM(n).SetLevelSet(u,Time,Level,t);     
        GFM(n).SetBndry(u,BF(u),Time,Level,t);
      }
    }
  }

  virtual void Initialize_(const double& dt_start) { 
    int* v = new int[NGFM()];
    for (register int n=0; n<NGFM(); n++) {
      v[n] = GFM(n).Boundary().GetVerbose();
      GFM(n).Boundary().SetVerbose(0);
    }
    base::Initialize_(dt_start);
    for (register int n=0; n<NGFM(); n++) 
      GFM(n).Boundary().SetVerbose(v[n]);
    delete [] v;
  }

  virtual void Output() {
    if (CurrentTime(base::GH(),0) == base::LastOutputTime()) return;
    if (base::_FileOutput) {
      START_WATCH 
        for (register int n=0; n<NGFM(); n++) 
          if (GFM(n).LevelSet().PlotPhi()) {
            double t = GetPhysicalTime(base::U(),CurrentTime(base::GH(),0),0);
            char name[DAGHBktGFNameWidth+16];
            std::sprintf(name,"phi_%d",n+1); 
            for (int l=0; l<=FineLevel(base::GH()); l++) {
              int Time = CurrentTime(base::GH(),l); 
              base::FileOutput_().WritePlain(GFM(n).Phi(), name, Time, l, t);  
            }
          }
      END_WATCH_OUPUT  
    }
    base::Output();
  }

  virtual bool Restart_(const char* CheckpointFile) {
    bool ret = base::Restart_(CheckpointFile);
    if (ret) {
      START_WATCH
        for (register int l=0; l<=FineLevel(base::GH()); l++) {
          int Time = CurrentTime(base::GH(),l);
          for (register int n=0; n<NGFM(); n++) {
            SetPhysicalTime(GFM(n).Phi(),Time,l,base::t[l]);
            if (base::ErrorEstimation()) {
              SetPhysicalTime(GFM(n).Phish(),Time,l,base::t[l]);
            }
          }
          _MaxRecomposeLevel = DAGHNull;
          SetRecomposeBndry(base::U(),Time,l,base::t[l]);
          if (base::ErrorEstimation()) 
            SetRecomposeBndry(base::Ush(),Time,l,base::t[l]);
        }
      END_WATCH_INITIALIZATION
    }
    return ret;
  }

  virtual void Checkpointing_(const char* CheckpointFile) {
    for (register int n=0; n<NGFM(); n++) {
      SetCheckpointFlag(GFM(n).Phi(), DAGHFalse);
      if (base::ErrorEstimation()) 
        SetCheckpointFlag(GFM(n).Phish(), DAGHFalse);
    }

    base::Checkpointing_(CheckpointFile);   

    for (register int n=0; n<NGFM(); n++) {
      SetCheckpointFlag(GFM(n).Phi(), DAGHTrue);
      if (base::ErrorEstimation()) 
        SetCheckpointFlag(GFM(n).Phish(), DAGHTrue);
    }
  }


  void AddGFM(gfm_type* gfm) {
    if (gfm) {
      gfm_type** _GFM_new = new gfm_type*[_nGFM+1];
      if (_GFM) {
        for (register int n=0; n<_nGFM; n++)
          _GFM_new[n] = _GFM[n];
        delete [] _GFM;
      }
      _GFM = _GFM_new;
      _GFM[_nGFM] = gfm;
      _nGFM++;
    }
  }

  void EliminateGFM(gfm_type* gfm) {
    if (gfm) {
      int c=0;
      for (register int n=0; n<_nGFM; n++)
        if (_GFM[n] == gfm) c++;
      if (c>0) {
        gfm_type** _GFM_new = (gfm_type**) 0;
        if (_nGFM-c>0) {
          _GFM_new = new gfm_type*[_nGFM-c];
          for (register int n=0, m=0; n<_nGFM; n++)
            if (_GFM[n] != gfm) 
              _GFM_new[m++] = _GFM[n];
        }
        delete [] _GFM;
        _GFM = _GFM_new;
        _nGFM = _nGFM-c;
      }
    }
  }

  void DeleteGFM(gfm_type* gfm) {
    EliminateGFM(gfm);
    if (gfm) {
      delete gfm->GetBoundaryP();
      delete gfm->GetLevelSetP();
      delete gfm;
    }
  }

  inline bool_grid_fct_type* BFP() { return _bf; }
  inline bool_grid_fct_type& BF() { return *_bf; }
  inline const bool_grid_fct_type& BF() const { return *_bf; }
  inline bool_grid_fct_type& BFsh() { return *_bf_sh; }
  inline const bool_grid_fct_type& BFsh() const{ return *_bf_sh; }
  inline bool_grid_fct_type& BF(vec_grid_fct_type& u) 
  { return (&u==base::_u_sh ? *_bf_sh : *_bf); }
  inline const bool_grid_fct_type& BF(vec_grid_fct_type& u) const 
  { return (&u==base::_u_sh ? *_bf_sh : *_bf); }

  inline gfm_type* GFMP(const int n) { assert (n>=0 && n<_nGFM); return _GFM[n]; }
  inline gfm_type& GFM(const int n) { assert (n>=0 && n<_nGFM); return *(_GFM[n]); }
  inline const gfm_type& GFM(const int n) const { assert (n>=0 && n<_nGFM); return *(_GFM[n]); }

  inline const int& NGFM() const { return _nGFM; }

protected:
  bool_grid_fct_type *_bf, *_bf_sh;
  gfm_recompose_functor_type* _GFMRecomposeFunc;
  gfm_type** _GFM;
  int _nGFM;
  int _MaxRecomposeLevel, _RecoverExterior;
  char BFName[DAGHBktGFNameWidth], BFNamesh[DAGHBktGFNameWidth];
};


#endif

---