ParticleAMRNS.H

Go to the documentation of this file.

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

// ParticleAMRNS.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, Fri, Jan 14, 2000

#ifndef _ParticleAMRNS_H_
#define _ParticleAMRNS_H_

#include "AMRLevel.H"
#include "LevelFluxRegister.H"
#include "CoarseAverage.H"
#include "PiecewiseLinearFillPatch.H"
#include "BinFab.H"
#include "DragParticle.H"
#include "ParticleProjector.H"
#include "CCProjector.H"
#include "LevelHelmholtzSolver.H"
#include "PoissonOp.H"
#include "VelocityPoissonOp.H"
#include "PhysBCUtil.H"


#ifdef NEW_ADVECTION
#include "PatchAdvection.H"
#include "PhysIBC.H"
#endif


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



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

  zMarker,
#endif

  NUM_SCAL_TYPES};




class ParticleAMRNS: public AMRLevel
{
public:
  
  ParticleAMRNS();
  
  ParticleAMRNS(AMRLevel* a_coarser_level_ptr, 
           const Box& a_prob_domain,
           int a_level, int a_ref_ratio);

  ParticleAMRNS(AMRLevel* a_coarser_level_ptr, 
           const ProblemDomain& a_prob_domain,
           int a_level, int a_ref_ratio);

  virtual ~ParticleAMRNS();


  virtual AMRLevel* makeAMRLevel(AMRLevel* a_coarser_level_ptr,
                                 const Box& a_problem_domain,
                                 int a_level, int a_ref_ratio) const;

  virtual AMRLevel* makeAMRLevel(AMRLevel* a_coarser_level_ptr,
                                 const ProblemDomain& a_problem_domain,
                                 int a_level, int a_ref_ratio) const;

  virtual void define(AMRLevel* a_coarse_level_ptr, 
                      const Box& a_problem_domain,
                      int a_level, int a_ref_ratio);

  virtual void define(AMRLevel* a_coarse_level_ptr, 
                      const ProblemDomain& a_problem_domain,
                      int a_level, int a_ref_ratio);

  virtual Real advance();

  virtual void postTimeStep();

  virtual void tagCells(IntVectSet& a_tags);

  virtual void tagCellsInit(IntVectSet& a_tags);

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



  virtual void postRegrid(int a_lBase);


  virtual void initialGrid(const Vector<Box>& a_new_grids);

  virtual void initialData();

  virtual void postInitialize();

#ifdef HDF5

  virtual void writeCheckpointHeader(HDF5Handle& a_handle) const;

  virtual void writeCheckpointLevel(HDF5Handle& a_handle) const;

  virtual void readCheckpointHeader(HDF5Handle& a_handle);

  virtual void readCheckpointLevel(HDF5Handle& a_handle);

  virtual void writePlotHeader(HDF5Handle& a_handle) const;

  virtual void writePlotLevel(HDF5Handle& a_handle) const;

  virtual void writeParticleData(HDF5Handle& a_handle) const;  

  virtual int writeParticlesHDF5(HDF5Handle& a_handle,
                                 const List<DragParticle>& a_particles) const;
#endif

  virtual Real computeDt();

  virtual Real computeInitialDt();

  virtual void CFL(Real a_cfl);

  virtual void particleCFL(Real a_particle_cfl);


  virtual void refinementThreshold(Real a_refine_threshold);

  void limitSolverCoarsening(bool a_limitSolverCoarsening);

  bool finestLevel() const;
  
  bool isEmpty() const;
  
  // compute vorticity
  void computeVorticity(LevelData<FArrayBox>& a_vorticity) const;

  // compute kinetic energy (only on interior cells)
  void computeKineticEnergy(LevelData<FArrayBox>& a_energy) const;

  // returns kinetic energy of all the particles in valid cells 
  Real particleKineticEnergy() const;

  // dumps ASCII list of particles
  void dumpParticlesASCII(std::ostream& os ) const;


  Real Dx() const;
  
  Real Nu() const;

  ParticleAMRNS* finerNSPtr() const;

  ParticleAMRNS* crseNSPtr() const;

  LevelData<FArrayBox>& newVel();

  const LevelData<FArrayBox>& newVel() const;

  LevelData<FArrayBox>& oldVel();

  const LevelData<FArrayBox>& oldVel() const;

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

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

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

  int numParticles() const;

  LevelData<FArrayBox>& newLambda();

  const LevelData<FArrayBox>& newLambda() const;

  LevelData<FArrayBox>& oldLambda();

  const LevelData<FArrayBox>& oldLambda() const;

  void lambda(LevelData<FArrayBox>& a_lambda, Real a_time) const;
  
  LevelData<FArrayBox>& newScal(const int a_comp);

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

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

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

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

  
  void setPhysBC(const PhysBCUtil& a_BC);

  PhysBCUtil* getPhysBCPtr() const;


protected:
  DisjointBoxLayout loadBalance(const Vector<Box>& a_grids);

  void levelSetup (const DisjointBoxLayout& level_domain);


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

  void readParameters();

  void finestLevel(bool a_finest_level);

  void swapOldAndNewStates();

  void resetStates(const Real a_time);

  void computeAdvectionVelocities(LevelData<FluxBox>& a_adv_vel,
                                  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.
  */
  void predictVelocities(LevelData<FArrayBox>& a_uDelU,
                         LevelData<FluxBox>& a_advVel,
                         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(LevelData<FArrayBox>& a_uStar,
                    LevelData<FArrayBox>& a_dragForce);



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


  void updateParticleForces(LevelData<BinFab<DragParticle> >& a_particles);

  void interpVelToParticles(LevelData<BinFab<DragParticle> >& a_particles,
                            const LevelData<FArrayBox>& a_velocity);
  
  void advectScalar(LevelData<FArrayBox>& a_new_scal,
                    LevelData<FArrayBox>& a_old_scal,
                    LevelData<FluxBox>& a_adv_vel,
                    LevelFluxRegister* a_crseFluxRegPtr,
                    LevelFluxRegister& a_fineFluxReg,
#ifdef NEW_ADVECTION
                    PhysIBC* a_scalIBC,
#endif
                    Real a_dt);

  void advectDiffuseScalar(LevelData<FArrayBox>& a_new_scal,
                           LevelData<FArrayBox>& a_old_scal,
                           LevelData<FluxBox>& a_adv_vel,
                           const Real& a_diffusive_coeff,
                           const LevelData<FArrayBox>* a_crse_scal_old,
                           const LevelData<FArrayBox>* a_crse_scal_new,
                           const Real a_old_crse_time,
                           const Real a_new_crse_time,
                           LevelFluxRegister* a_crseFluxRegPtr,
                           LevelFluxRegister& a_fineFluxReg,
                           DomainGhostBC& a_scalPhysBC,
#ifdef NEW_ADVECTION
                    PhysIBC* a_scalIBC,
#endif
                           Real a_dt);

  void initializeGlobalPressure();

  void initializeLevelPressure(Real a_currentTime, Real a_dtInit);

  void smoothVelocityField(int lbase);

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

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

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

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

  void fillLambda(LevelData<FArrayBox>& a_lambda, Real a_time) const;

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


  virtual int numPlotComps() const;

  virtual void getPlotData(LevelData<FArrayBox>& a_plot_data) const;

  void  defineRegridAMRSolver(AMRSolver& a_solver, 
                              const Vector<DisjointBoxLayout>& a_grids,
                              const Vector<ProblemDomain>& a_domains,
                              const Vector<Real>& a_amrDx,
                              const Vector<int>& a_refRatios,
                              const int& a_lBase);
  




protected:
  LevelData<FArrayBox>* m_vel_old_ptr;

  LevelData<FArrayBox>* m_vel_new_ptr;

#ifdef DEBUG

  LevelData<FArrayBox>* m_vel_save_ptr;  

  Real m_saved_time;

#endif

  LevelData<BinFab<DragParticle> > m_particles;

  LevelData<FArrayBox>* m_lambda_old_ptr;
  
  LevelData<FArrayBox>* m_lambda_new_ptr;

#ifdef DEBUG

  LevelData<FArrayBox>* m_lambda_save_ptr;  
#endif

  Vector<LevelData<FArrayBox>*> m_scal_old;
  
  Vector<LevelData<FArrayBox>*> m_scal_new;

#ifdef DEBUG

  Vector<LevelData<FArrayBox>*> m_scal_save;
#endif

  static const int s_num_vel_comps = SpaceDim;

  static int s_num_scal_comps;

  static const int s_max_scal_comps=SpaceDim;

  static const char* s_vel_names[s_num_vel_comps];

  static const char* s_scal_names[s_max_scal_comps];

  static Vector<Real> s_scal_coeffs;

  static Vector<Real> s_domLength;
  
  static bool s_ppInit;

  static Real s_init_shrink;

  static Real s_max_dt;

  static Real s_prescribedDt;

  static bool s_project_initial_vel;

  static bool s_initialize_pressures;

  static bool s_smooth_after_regrid;

  static bool s_dump_smoothing_plots;

  static Real s_regrid_smoothing_coeff;

  static int s_num_init_passes;

  static bool s_reflux_momentum;

  static bool s_reflux_normal_momentum;

  static bool s_set_bogus_values;

  static Real s_bogus_value;

  static bool s_tag_vorticity;

  static Real s_vort_factor;

  static int s_tags_grow;

  static bool s_implicit_reflux;

  static bool s_applyFreestreamCorrection;

  static bool s_reflux_scal;

  static bool s_implicit_scal_reflux;

  static int s_viscous_solver_type;

  static Real s_viscous_solver_tol;
  
  static int s_viscous_num_smooth_up;

  static int s_viscous_num_smooth_down;

  static Real s_nu;

  static bool s_do_particles;

  static bool s_read_particles_from_file;

  static string s_particle_file;

  static Real s_particle_drag_coeff;

  static RealVect s_particle_body_force;

  static Real s_delta_epsilon;

  static bool s_use_image_particles;

  static int s_order_vel_interp;

  static bool s_specifyInitialGrids;

  static string s_initialGridFile;

  static bool s_init_vel_from_vorticity;
  
  static Real s_backgroundVel;


  static bool s_write_particles;
  
  static bool s_write_divergence;

  static bool s_write_lambda;

  static bool s_write_time_derivatives;

  static bool s_write_vorticity;

  static bool s_write_scalars;

  static bool s_write_dScalar_dt;

  static bool s_write_strains;

  static bool s_write_grad_eLambda;  

  static bool s_write_error;

  //static bool s_write_curl_error;

  static bool s_write_proc_ids;

  static bool s_compute_scal_err;

  static bool s_write_initial_solve_rhs;

  static bool s_write_grids;


  Vector<Tuple<Real, SpaceDim> > m_dataPoints;

  Real m_dx;

  Real m_cfl;

  Real m_particleCFL;

  bool m_limitSolverCoarsening;


  Real m_dt_save;

  int m_level_steps;

  CoarseAverage m_coarse_average;

  CoarseAverage m_coarse_average_lambda;

#if 0

  CoarseAverage m_coarse_average_scal;
#endif

  CCProjector m_projection;
  
  ParticleProjector m_particle_projector;

  LevelFluxRegister m_flux_register;

  LevelFluxRegister m_lambda_flux_reg;

  Vector<LevelFluxRegister*> m_scal_fluxreg_ptrs;

  QuadCFInterp m_velCFInterp; 

  VelocityPoissonOp m_velocityPoissonOp;

  PoissonOp m_scalarsPoissonOp;

  LevelHelmholtzSolver m_viscous_solver;

#ifdef NEW_ADVECTION
  PatchAdvection m_patchGod;
#endif

  Real m_refine_threshold;

  bool m_finest_level;

  bool m_is_empty;

  bool m_regrid_smoothing_done;

  PhysBCUtil* m_physBCPtr;



};

#endif

---