---

PAmrNS.H

Go to the documentation of this file.

/* _______              __
  / ___/ /  ___  __ _  / /  ___
 / /__/ _ \/ _ \/  ' \/ _ \/ _ \
 \___/_//_/\___/_/_/_/_.__/\___/ 
*/

// PAmrNS.H

//
// This software is copyright (C) by the Lawrence Berkeley
// National Laboratory.  Permission is granted to reproduce
// this software for non-commercial purposes provided that
// this notice is left intact.
// 
// It is acknowledged that the U.S. Government has rights to
// this software under Contract DE-AC03-765F00098 between
// the U.S.  Department of Energy and the University of
// California.
//
// This software is provided as a professional and academic
// contribution for joint exchange. Thus it is experimental,
// is provided ``as is'', with no warranties of any kind
// whatsoever, no support, no promise of updates, or printed
// documentation. By using this software, you acknowledge
// that the Lawrence Berkeley National Laboratory and
// Regents of the University of California shall have no
// liability with respect to the infringement of other
// copyrights by any part of this software.
//
// Dan Martin, Sept 13, 2004

#ifndef _PAmrNS_H_
#define _PAmrNS_H_

#include "AMRLevel.H"
#include "CoarseAverage.H"
#include "PiecewiseLinearFillPatch.H"
#include "BinFab.H"
#include "DragParticle.H"
#include "AMRParticleProjector.H"
#include "AmrProjector.H"
#include "PoissonOp.H"
#include "VelocityPoissonOp.H"
#include "PhysBCUtil.H"


#include "PatchAdvection.H"
#include "OldPhysIBC.H"

#define ADVECT_GROW 4

#ifdef CH_USE_DOUBLE
#define TIME_EPS 1.0e-10
#else
#define TIME_EPS 1.0e-5
#endif

// small parameter used in defining TGA constants for viscous solves
#define TGA_EPS 1.0e-8


enum viscousSolverTypes {
  backwardEuler = 0,
  CrankNicolson,
  TGA, 
  NUM_SOLVER_TYPES};


enum particleUpdateTypes {
  predictorCorrector = 0,
  forwardEuler,
  betterPredictorCorrector,
  betterForwardEuler, 
  NUM_UPDATE_TYPES};


enum scalTypes {
  xMarker = 0,
  yMarker,
#if CH_SPACEDIM >= 3

  zMarker,
#endif

  NUM_SCAL_TYPES};


// number to use to compare real numbers
#define REALSMALL 1.0e-10



class PAmrNS
{
public:
  
  PAmrNS();
  

  PAmrNS(const ProblemDomain& a_vectDomain,
         const Vector<int>& a_vectRefRatio,
         const RealVect& a_domainLength,
         int a_maxLevel);


  ~PAmrNS();

  

  void define(const ProblemDomain& a_baseDomain,
              const Vector<int>& a_vectRefRatio,
              const RealVect& a_domainLength, 
              int a_maxLevel);
         
  void setupForFixedHierarchyRun(Vector<Vector<Box> >& a_vectBoxes);

  void setupForAmrRun();

  void setupForRestart(const string& a_restartFile);

  void setupForRestart(const HDF5Handle& a_restartFile);

  Real doRun(Real a_maxTime, int a_maxStep);

  void conclude();

#ifdef HDF5

  int writeCheckpointFile();
  
  int writeCheckpointFile(const string& a_fname);

  int writeCheckpointFileHDF5(HDF5Handle& a_handle);
  
  int readCheckpointFile(const string& a_fname);

  int readCheckpointFile(HDF5Handle& a_handle);

  int writePlotFile();

  int writePlotFileName(const string& a_fname);

  int writePlotFileHDF5(HDF5Handle& a_handle);
#endif


  void checkpointInterval(int a_checkpointInterval);
  
  void plotInterval(int a_plotInterval);

  void checkpointPrefix(const string& a_prefix);

  void plotPrefix(const string& a_prefix);
  
  void maxGridSize(int a_maxGridSize);

  void maxBaseGridSize(int a_maxBaseGridSize);

  void fillRatio(Real a_fillRatio);

  void blockFactor(int a_blockFactor);

  void regridInterval(int a_regridInterval);

  void maxDtGrow(Real& a_maxDtGrow);

  void fixedDt(Real a_dt); 
    
  void CFL(Real a_cfl);
  
  void particleCFL(Real a_particle_cfl);
  
  
  void refinementThreshold(Real a_refine_threshold);
  
  void limitSolverCoarsening(bool a_limitSolverCoarsening);

  void verbosity(int a_verbosity); 

  
  
  int finestLevel() const {return m_finestLevel;};
  
  int maxLevel() const {return m_maxLevel;};
  
  bool isEmpty(int a_level) const {return (a_level > m_finestLevel);};
  

  Real Dx(int a_level) const {return m_vectDx[a_level];};
  
  Real Nu() const {return m_nu;};

  void setPhysBC(const PhysBCUtil& a_BC);


  bool isDefined() const;


protected:


  void advance(Real a_dt);
  
  void postTimeStep();
  
  void tagCells(Vector<IntVectSet>& a_tags);
  
  void tagCellsLevel(IntVectSet& a_tags, int a_level);
  
  void tagCellsInit(Vector<IntVectSet>& a_tags);
  

  void regrid();

  
  void postRegrid();
  
  

  void initialGrid(const Vector<Vector<Box> >& a_new_grids,
                   int a_baseLevel,
                   bool a_clearStorage = false);



  void initialGrid(const Vector<DisjointBoxLayout>& a_new_grids,
                   const Vector<DisjointBoxLayout>& a_new_particle_grids,
                   int a_baseLevel,
                   bool a_clearStorage = false);
  
  void initialData();
  
  void postInitialize();
  
#ifdef HDF5

  void writeParticleData(HDF5Handle& a_handle) const;  
  
  int writeParticlesHDF5(HDF5Handle& a_handle,
                         const List<DragParticle>& a_particles) const;
#endif
  
  Real computeDt();
  
  Real computeInitialDt();
  
  void computeVorticity(LevelData<FArrayBox>& a_vorticity, int a_level) const;

  void computeKineticEnergy(LevelData<FArrayBox>& a_energy,
                            int a_level) const;

  Real particleKineticEnergy() const;

  void dumpParticlesASCII(std::ostream& os ) const;


  void dumpLoadBalanceInfo() const;

  LevelData<FArrayBox>& newVel(int a_level);

  const LevelData<FArrayBox>& newVel(int a_level) const;

  LevelData<FArrayBox>& oldVel(int a_level);

  const LevelData<FArrayBox>& oldVel(int a_level) const;

  void velocity(LevelData<FArrayBox>& a_vel, Real a_time, int a_level) const;

  LevelData<BinFab<DragParticle> >& particles(int a_level);

  const LevelData<BinFab<DragParticle> >& particles(int a_level) const;

  int numParticles(int a_level) const;
  
  LevelData<FArrayBox>& newScal(const int a_comp, int a_level);

  const LevelData<FArrayBox>& newScal(const int a_comp, int a_level) const;

  LevelData<FArrayBox>& oldScal(const int a_comp, int a_level);

  const LevelData<FArrayBox>& oldScal(const int a_comp, int a_level) const;

  void scalars(LevelData<FArrayBox>& a_scalars, const Real a_time, 
               const int a_comp, int a_level) const;

  
  PhysBCUtil* getPhysBCPtr() const;

  

protected:
  void loadBalance(Vector<DisjointBoxLayout>& a_vectGrids,
                   Vector<DisjointBoxLayout>& a_vectParticleGrids,
                   const Vector<Vector<Box> >& a_vectBoxes,
                   int a_lbase); 

  void setup ();

  void regrid(const Vector<Vector<Box> >& a_new_grids);
  


  void defineParticles(List<DragParticle>& a_particles) const;

  void readParameters();

  void finestLevel(int a_finest_level);

  void swapOldAndNewStates();

  void resetStates(const Real a_time);


  void computeAdvectionVelocities(Vector<LevelData<FluxBox>*> & a_adv_vel,
                                  Vector<LevelData<FArrayBox>* >& a_phi,
                                  Vector<LevelData<FArrayBox>* >& a_dragForce);


  /*@ManMemo: do predictor -- returns approx. to (U dot del)U.
    Also updates flux registers with momentum fluxes and
    part of viscous fluxes which contains the old-time velocity.
    a_phi is the correction from the MAC projection
  */
  void predictVelocities(Vector<LevelData<FArrayBox>* >& a_uDelU,
                         Vector<LevelData<FluxBox>* >& a_advVel,
                         Vector<LevelData<FArrayBox>* >& a_phi,
                         Vector<LevelData<FArrayBox>* >& a_dragForce);
  
  /*@ManMemo: This function computes uStar -- the approximation
    of the new-time velocity which will be projected.  The
    advection term (u-dot-del-U) is passed in in a_uStar, which 
    is modified in place to contain uStar. Note that both viscous
    and inviscid cases, this returns the approximation to 
    u^{n+1} + dt*gradPi.  a_dragForce is the old-time drag force
  */
  void computeUStar(Vector<LevelData<FArrayBox>* >& a_uStar,
                    Vector<LevelData<FArrayBox>* >& a_dragForce,
                    Vector<LevelData<FArrayBox>* >& a_Pi);



  void updateParticlePositions(Vector<LevelData<FArrayBox>*>& a_newDrag,
                               Vector<LevelData<BinFab<DragParticle> >* >& a_particles,
                               const Vector<LevelData<FArrayBox>*>& a_fluidVel,
                               bool a_completeUpdate);


  void updateParticleForces(Vector<LevelData<BinFab<DragParticle> >* >& a_particles);
  
  void interpVelToParticles(Vector<LevelData<BinFab<DragParticle> >* >& a_particles,
                            const Vector<LevelData<FArrayBox>* >& a_velocity);
  
  void advectScalar(Vector<LevelData<FArrayBox>* >& a_new_scal,
                    Vector<LevelData<FArrayBox>* >& a_old_scal,
                    Vector<LevelData<FluxBox>* >& a_adv_vel,
                    OldPhysIBC* a_scalIBC,
                    Real a_dt);

  void advectDiffuseScalar(Vector<LevelData<FArrayBox>* >& a_new_scal,
                           Vector<LevelData<FArrayBox>* >& a_old_scal,
                           Vector<LevelData<FluxBox>* >& a_adv_vel,
                           const Real& a_diffusive_coeff,
                           DomainGhostBC& a_scalPhysBC,
                           OldPhysIBC* a_scalIBC,
                           Real a_dt);

  void initializePressure();

  void computeLapVel(Vector<LevelData<FArrayBox>* >& a_lapVel,
                     Vector<LevelData<FArrayBox>* >& a_vel);

  void computeLapScal(Vector<LevelData<FArrayBox>* >& a_lapScal,
                      Vector<LevelData<FArrayBox>* >& a_scal,
                      DomainGhostBC& a_physBC);


  void fillVelocity(Vector<LevelData<FArrayBox>* >& a_vel, Real a_time);

  void fillVelocity(Vector<LevelData<FArrayBox>* >& a_vel, Real a_time,  
                    int vel_comp, int dest_comp, int num_comp);

  void fillScalars(Vector<LevelData<FArrayBox>* >& a_scal, Real a_time,
                   const int a_comp) const;


  int numPlotComps() const;
  
  void getPlotDataNames(Vector<string>& a_vectNames);
  
  void getPlotData(Vector<LevelData<FArrayBox>* >& a_plot_data);

  void fillCheckpointData(Vector<LevelData<FArrayBox>* >& a_vectData,
                          Vector<string>& a_checkDataNames);

  void setDefaultValues();


protected:
  Vector<int> m_vectRefRatio;
  
  Vector<ProblemDomain> m_vectDomain;

  Vector<Real> m_vectDx;

  Vector<DisjointBoxLayout> m_vectGrids;

  Vector<DisjointBoxLayout> m_vectParticleGrids;

  Vector<Vector<Box> > m_vectBoxes;


  int m_maxLevel;

  int m_finestLevel;

  int m_regrid_interval;

  int m_block_factor;

  Real m_fill_ratio;
  
  int m_nesting_radius;
  
  int m_max_box_size;

  int m_max_base_box_size;

  Real m_time;

  int m_cur_step;

  Real m_dt;

  Vector<LevelData<FArrayBox>* > m_vectOldVel;

  Vector<LevelData<FArrayBox>* > m_vectNewVel;

  Vector<LevelData<FArrayBox>* > m_vectPi;

  Vector<LevelData<BinFab<DragParticle> >* > m_vectParticles;


  Vector<Vector<LevelData<FArrayBox>* > > m_vectOldScal;
  
  Vector<Vector<LevelData<FArrayBox>* > > m_vectNewScal;


  int m_num_vel_comps;
  
  int m_num_scal_comps;

  int m_max_scal_comps;

  static const char* s_vel_names[CH_SPACEDIM];

  static const char* s_scal_names[NUM_SCAL_TYPES];
  
  Vector<Real> m_scal_coeffs;
  
  RealVect m_domLength;
    
  Real m_init_shrink;

  Real m_max_dt;

  Real m_min_dt;
  
  Real m_prescribedDt;

  Real m_max_dt_grow;

  int m_verbosity;

  bool m_project_initial_vel;

  bool m_update_velocity;
  
  bool m_initialize_pressures;
  

  int m_num_init_passes;
  
  bool m_set_bogus_values;

  Real m_bogus_value;

  bool m_tag_vorticity;

  Real m_vort_factor;

  bool m_tag_particles;
  
  int m_tags_grow;

  int m_viscous_solver_type;

  Real m_viscous_solver_tol;
  
  int m_viscous_num_smooth_up;

  int m_viscous_num_smooth_down;

  bool m_limitSolverCoarsening; 

  Real m_nu;

  bool m_doRestart;

  bool m_do_particles;

  bool m_read_particles_from_file;

  string m_particle_file;
  
  Real m_particle_drag_coeff;

  RealVect m_particle_body_force;

  Real m_particle_epsilon;

  Real m_spreading_radius;

  Real m_correction_radius;


  int m_particle_grids_buffer;

  bool m_use_image_particles;

  int m_order_vel_interp;

  bool m_limit_vel_interp;

  int m_particle_update_type;
  
  bool m_specifyInitialGrids;

  string m_initialGridFile;

  bool m_init_vel_from_vorticity;
  
  Real m_backgroundVel;

  
  bool m_write_particles;
  
  bool m_write_divergence;

  bool m_write_time_derivatives;
  
  bool m_write_vorticity;

  bool m_write_scalars;

  bool m_write_dScalar_dt;
  
  bool m_write_error;

  // bool m_write_curl_error;

  bool m_write_proc_ids;

  bool m_compute_scal_err;

  bool m_write_initial_solve_rhs;
  
  bool m_write_grids;

  Real m_cfl;

  Real m_particleCFL;


  Vector<CoarseAverage*> m_coarse_average;
  
  Vector<CoarseAverage*> m_coarse_average_scal;

  AmrProjector m_projection;
  
  AMRParticleProjector m_particle_projector;
 
  int m_particleProjCoarsening;

  Vector<QuadCFInterp*> m_velCFInterp_ptrs;     

  Vector<PoissonOp*> m_scalarsPoissonOp_ptrs;
  
  Vector<PatchAdvection*> m_patchGod;

  Real m_refine_threshold;
  
  PhysBCUtil* m_physBCPtr;


  int m_plot_interval;
  
  int m_check_interval;


  int m_lastcheck_step;


  int m_restart_step;      

  Vector<int> m_cell_updates;

  bool m_ppInit;

  bool m_isDefined;
  

  std::string m_plotfile_prefix;
  std::string m_checkpointfile_prefix;

};

#endif

---