Chombo + EB + MF  3.2
EBPoissonOp.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _EBPOISSONOP_H__
12 #define _EBPOISSONOP_H__
13 
14 #include "REAL.H"
15 #include "Box.H"
16 #include "FArrayBox.H"
17 #include "Vector.H"
18 #include <map>
19 #include "RefCountedPtr.H"
20 
21 #include "AMRMultiGrid.H"
22 
23 #include "EBIndexSpace.H"
24 #include "EBCellFAB.H"
25 #include "EBCellFactory.H"
26 
27 #include "EBStencil.H"
28 
29 #include "EBLevelDataOps.H"
30 #include "BaseEBBC.H"
31 #include "BaseDomainBC.H"
32 #include "CFIVS.H"
33 #include "EBFluxRegister.H"
34 #include "EBMGAverage.H"
35 #include "EBMGInterp.H"
36 #include "EBCoarsen.H"
37 #include "PolyGeom.H"
38 #include "EBQuadCFInterp.H"
39 #include "EBLevelGrid.H"
40 #include "NamespaceHeader.H"
41 
42 
43 // Weight for weighted Jacobi. This kind of ought to be a member
44 // variable that could be set to a different value at run time. For
45 // now, this gives the best convergence.
46 #define EBPOISSONOP_JACOBI_OMEGA 0.67
47 
48 
49 #if CH_SPACEDIM==2
50 #define EBPO_NUMSTEN 4
51 #elif CH_SPACEDIM==3
52 #define EBPO_NUMSTEN 8
53 #else
55 {
56  THIS_WILL_ONLY_COMPILE_WHEN_CH_SPACEDIM_IS_2_OR_3;
57 }
58 #endif
59 ///
60 /**
61  Operator to solve (alpha + beta lapl)phi = rhs. This follows the AMRLevelOp interface.
62 */
63 class EBPoissonOp: public MGLevelOp<LevelData<EBCellFAB> >
64 {
65 public:
66  ///
67  ~EBPoissonOp();
68 
69  ///
70  EBPoissonOp();
71 
72  ///
73  /**
74  If you are approaching this operator from this interface, consider backing away and using
75  EBPoissonOpFactory to generate these objects. Really.
76  a_eblg, : grid at this level \\
77  a_eblgCoarMG, : grid at intermediate multigrid level \\
78  a_domainBC, : domain boundary conditions at this level \\
79  a_ebBC: eb boundary conditions at this level \\
80  a_dx: grid spacing at this level \\
81  a_origin: offset to lowest corner of the domain \\
82  a_hasCoarserMG: true if there is a coarser MultiGrid level. \\
83  a_preCondIters: number of iterations to do for pre-conditioning \\
84  a_relaxType: 0 means point Jacobi, 1 is Gauss-Seidel. \\
85  a_orderEB: 0 to not do flux interpolation at cut faces. \\
86  a_alpha: coefficent of identity \\
87  a_beta: coefficient of laplacian.\\
88  a_ghostCellsPhi: Number of ghost cells in phi, correction (typically one)\\
89  a_ghostCellsRhs: Number of ghost cells in RHS, residual, lphi (typically zero)\\
90  Ghost cell arguments are there for caching reasons. Once you set them, an error is thrown if
91  you send in data that does not match.
92  */
93  EBPoissonOp( const EBLevelGrid & a_eblg,
94  const EBLevelGrid & a_eblgCoarMG,
95  const RefCountedPtr<BaseDomainBC>& a_domainBC,
96  const RefCountedPtr<BaseEBBC>& a_ebBC,
97  const RealVect& a_dx,
98  const RealVect& a_origin,
99  const bool& a_hasMGObjects,
100  const int& a_numPreCondIters,
101  const int& a_relaxType,
102  const int& a_orderEB,
103  const Real& a_alpha,
104  const Real& a_beta,
105  const IntVect& a_ghostCellsPhi,
106  const IntVect& a_ghostCellsRHS);
107 
108  //MGOp operations. no finer or coarser
109 
110  ///
111  /**
112  */
113  virtual void residual(LevelData<EBCellFAB>& a_residual,
114  const LevelData<EBCellFAB>& a_phi,
115  const LevelData<EBCellFAB>& a_rhs,
116  bool a_homogeneousPhysBC=false);
117 
118  ///
119  /**
120  */
121  virtual void preCond(LevelData<EBCellFAB>& a_opPhi,
122  const LevelData<EBCellFAB>& a_phi);
123 
124 
125  virtual void
127  const LevelData<EBCellFAB>& a_rhs,
128  const IntVect& a_color,
129  const Real& a_weight,
130  const bool& a_homogeneousPhysBC);
131 
132  virtual void
134  const LevelData<EBCellFAB>& a_rhs,
135  const IntVect& a_color,
136  const bool& a_homogeneousPhysBC,
137  int icolor);
138 
139 
140  ///
141  /**
142  this is the linearop function.
143  */
144  virtual void applyOp(LevelData<EBCellFAB>& a_opPhi,
145  const LevelData<EBCellFAB>& a_phi,
146  bool a_homogeneousPhysBC,
147  DataIterator& a_dit,
148  bool do_exchange = true);
149 
150  virtual void applyOp(LevelData<EBCellFAB>& a_opPhi,
151  const LevelData<EBCellFAB>& a_phi,
152  bool a_homogeneousPhysBC)
153  {
154  DataIterator dit = a_opPhi.dataIterator();
155  applyOp(a_opPhi, a_phi, a_homogeneousPhysBC, dit, true);
156  }
157 
158  ///
159  /**
160  */
161  virtual void create(LevelData<EBCellFAB>& a_lhs,
162  const LevelData<EBCellFAB>& a_rhs);
163 
164  ///
165  /**
166  */
167  virtual void assign(LevelData<EBCellFAB>& a_lhs,
168  const LevelData<EBCellFAB>& a_rhs);
169 
170  ///
171  /**
172  */
173  virtual Real dotProduct(const LevelData<EBCellFAB>& a_1,
174  const LevelData<EBCellFAB>& a_2);
175 
176  ///
177  /**
178  */
179  virtual void incr(LevelData<EBCellFAB>& a_lhs,
180  const LevelData<EBCellFAB>& a_x,
181  Real a_scale);
182 
183  ///
184  /**
185  */
186  virtual void axby(LevelData<EBCellFAB>& a_lhs,
187  const LevelData<EBCellFAB>& a_x,
188  const LevelData<EBCellFAB>& a_y,
189  Real a_a,
190  Real a_b);
191 
192  ///
193  /**
194  */
195  virtual void scale(LevelData<EBCellFAB>& a_lhs,
196  const Real& a_scale);
197 
198  ///
199  /**
200  */
201  virtual Real norm(const LevelData<EBCellFAB>& a_rhs,
202  int a_ord);
203 
204  ///
205  /**
206  */
207  virtual void setToZero(LevelData<EBCellFAB>& a_lhs);
208 
209  ///
210  /**
211  */
212  virtual void setVal(LevelData<EBCellFAB>& a_lhs, const Real& a_value);
213 
214  ///
215  /**
216  */
217  virtual void createCoarser(LevelData<EBCellFAB>& a_coarse,
218  const LevelData<EBCellFAB>& a_fine,
219  bool a_ghosted);
220 
221  ///
222  /**
223  */
224  virtual void relax(LevelData<EBCellFAB>& a_e,
225  const LevelData<EBCellFAB>& a_residual,
226  int a_iterations);
227 
228  ///
229  /**
230  Calculate restricted residual:
231  a_resCoarse[2h] = I[h->2h] (a_rhsFine[h] - L[h](a_phiFine[h]))
232  */
233  virtual void restrictResidual(LevelData<EBCellFAB>& a_resCoarse,
234  LevelData<EBCellFAB>& a_phiFine,
235  const LevelData<EBCellFAB>& a_rhsFine);
236 
237  ///
238  /**
239  Correct the fine solution based on coarse correction:
240  a_phiThisLevel += I[2h->h] (a_correctCoarse)
241  */
242  virtual void prolongIncrement(LevelData<EBCellFAB>& a_phiThisLevel,
243  const LevelData<EBCellFAB>& a_correctCoarse);
244 
245 
246  static bool nextColor(IntVect& color,
247  const IntVect& limit);
248 
249  static int getIndex(const IntVect& a_color);
250 
251  static IntVect getColor(const int& a_icolor);
252 
254  const LevelData<EBCellFAB>& a_rhs);
255 
256  void colorGS(LevelData<EBCellFAB>& a_phi,
257  const LevelData<EBCellFAB>& a_rhs,
258  const IntVect& color, int icolor);
259 
260 
261 
262 protected:
263 
264  void defineStencils();
265 
269 
270 
273 
276 
278 
289 
290 
291  //stencils for operator evaluation
293  //stencils for operator evaluation on gauss-seidel colors
295  //stencils for lambda*rhs addition to opphi
297  //weights that get multiplied by alpha
299  //weights that get multiplied by beta
301 
302  //cache the irreg vofiterator
305 
308 
311 
313  //stuff below is only defined if m_hasMGObjects==true
319 
320  /// Get the coefficients (for all cells) to multiply the residual by
321  /// before adding it to the previous iterate for Jacobi iteration.
322  void getJacobiRelaxCoeff(LevelData<EBCellFAB>& a_relaxCoeff);
323 
324 private:
325 
326 
327  //internally useful but not for general consumption
328  //lots of hidden assumptions and all that
329 
330  void getOpVoFStencil(VoFStencil& a_stencil,
331  const EBISBox& a_ebisbox,
332  const VolIndex& a_vof);
333 
334  void getOpVoFStencil(VoFStencil& a_stencil,
335  const int& a_dir,
336  const Vector<VolIndex>& a_allMonotoneVoFs,
337  const EBISBox& a_ebisbox,
338  const VolIndex& a_vof,
339  const bool& a_lowOrder);
340 
341 
342  void getOpFaceStencil(VoFStencil& a_stencil,
343  const Vector<VolIndex>& a_allMonotoneVofs,
344  const EBISBox& a_ebisbox,
345  const VolIndex& a_vof,
346  int a_dir,
347  const Side::LoHiSide& a_side,
348  const FaceIndex& a_face,
349  const bool& a_lowOrder);
350 
351 
352  void levelJacobi(LevelData<EBCellFAB>& a_phi,
353  const LevelData<EBCellFAB>& a_rhs);
354 
355  ///
356  /** applyOp() on the regular cells. Only used in line relaxation stuff
357  opPhi comes in holding alpha*phi. this adds on beta*lapl phi*/
358  void applyOpRegular( int idir,
359  Box * a_loBox,
360  Box * a_hiBox,
361  int * a_hasLo,
362  int * a_hasHi,
363  Box & a_curOpPhiBox,
364  Box & a_curPhiBox,
365  int a_nComps,
366  BaseFab<Real> & a_curOpPhiFAB,
367  const BaseFab<Real> & a_curPhiFAB,
368  bool a_homogeneousPhysBC,
369  const DataIndex& a_dit,
370  const Real& a_beta);
371 
372  ///
373  /** applyOp() on the regular cells for all directions
374  opPhi comes in holding alpha*phi. this adds on beta*lapl phi*/
375  void applyOpRegularAllDirs(Box * a_loBox,
376  Box * a_hiBox,
377  int * a_hasLo,
378  int * a_hasHi,
379  Box & a_curOpPhiBox,
380  Box & a_curPhiBox,
381  int a_nComps,
382  BaseFab<Real> & a_curOpPhiFAB,
383  const BaseFab<Real> & a_curPhiFAB,
384  bool a_homogeneousPhysBC,
385  const DataIndex& a_dit,
386  const Real& a_beta);
387 
388  void applyDomainFlux(Box * a_loBox,
389  Box * a_hiBox,
390  int * a_hasLo,
391  int * a_hasHi,
392  Box & a_curPhiBox,
393  int a_nComps,
394  BaseFab<Real> & a_phiFAB,
395  bool a_homogeneousPhysBC,
396  const DataIndex& a_dit,
397  const Real& a_beta);
398 
399  ///
400  /**
401  This function computes: a_lhs = (1/diagonal)*a_rhs
402  Use this function to initialize the preconditioner
403  */
405  const LevelData<EBCellFAB>& a_rhs);
406 
407 private:
408  //copy constructor and operator= disallowed for all the usual reasons
409  EBPoissonOp(const EBPoissonOp& a_opin)
410  {
411  MayDay::Error("invalid operator");
412  }
413 
414  void operator=(const EBPoissonOp& a_opin)
415  {
416  MayDay::Error("invalid operator");
417  }
418 };
419 
420 
421 #include "NamespaceFooter.H"
422 #endif
virtual void GSColorAllRegular(LevelData< EBCellFAB > &a_phi, const LevelData< EBCellFAB > &a_rhs, const IntVect &a_color, const Real &a_weight, const bool &a_homogeneousPhysBC)
EBISLayout m_ebislCoarMG
Definition: EBPoissonOp.H:317
const IntVect m_ghostCellsPhi
Definition: EBPoissonOp.H:266
Real m_dxScale
Definition: EBPoissonOp.H:281
Real m_beta
Definition: EBPoissonOp.H:283
LayoutData< VoFIterator > m_vofItIrregColorDomLo[EBPO_NUMSTEN][SpaceDim]
Definition: EBPoissonOp.H:309
int m_numPreCondIters
Definition: EBPoissonOp.H:287
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:141
RealVect m_origin
Definition: EBPoissonOp.H:285
void operator=(const EBPoissonOp &a_opin)
Definition: EBPoissonOp.H:414
void applyDomainFlux(Box *a_loBox, Box *a_hiBox, int *a_hasLo, int *a_hasHi, Box &a_curPhiBox, int a_nComps, BaseFab< Real > &a_phiFAB, bool a_homogeneousPhysBC, const DataIndex &a_dit, const Real &a_beta)
void colorGS(LevelData< EBCellFAB > &a_phi, const LevelData< EBCellFAB > &a_rhs, const IntVect &color, int icolor)
virtual void incr(LevelData< EBCellFAB > &a_lhs, const LevelData< EBCellFAB > &a_x, Real a_scale)
const IntVect m_ghostCellsRHS
Definition: EBPoissonOp.H:267
virtual void preCond(LevelData< EBCellFAB > &a_opPhi, const LevelData< EBCellFAB > &a_phi)
Definition: FaceIndex.H:28
Data that maintains a one-to-one mapping of T to the boxes in a BoxLayout.
Definition: BoxLayout.H:26
LayoutData< RefCountedPtr< EBStencil > > m_colorEBStencil[EBPO_NUMSTEN]
Definition: EBPoissonOp.H:294
Definition: EBISBox.H:46
Definition: EBLevelGrid.H:30
DataIterator dataIterator() const
Definition: LayoutDataI.H:78
virtual void setToZero(LevelData< EBCellFAB > &a_lhs)
static bool nextColor(IntVect &color, const IntVect &limit)
void getOpFaceStencil(VoFStencil &a_stencil, const Vector< VolIndex > &a_allMonotoneVofs, const EBISBox &a_ebisbox, const VolIndex &a_vof, int a_dir, const Side::LoHiSide &a_side, const FaceIndex &a_face, const bool &a_lowOrder)
LayoutData< VoFIterator > m_vofItIrreg
Definition: EBPoissonOp.H:303
Definition: DataIterator.H:190
void levelMulticolorGS(LevelData< EBCellFAB > &a_phi, const LevelData< EBCellFAB > &a_rhs)
virtual void scale(LevelData< EBCellFAB > &a_lhs, const Real &a_scale)
virtual void applyOp(LevelData< EBCellFAB > &a_opPhi, const LevelData< EBCellFAB > &a_phi, bool a_homogeneousPhysBC)
Definition: EBPoissonOp.H:150
void defineStencils()
virtual void axby(LevelData< EBCellFAB > &a_lhs, const LevelData< EBCellFAB > &a_x, const LevelData< EBCellFAB > &a_y, Real a_a, Real a_b)
LayoutData< VoFIterator > m_vofItIrregDomHi[SpaceDim]
Definition: EBPoissonOp.H:307
LayoutData< VoFIterator > m_vofItIrregColor[EBPO_NUMSTEN]
Definition: EBPoissonOp.H:304
LayoutData< RefCountedPtr< EBStencil > > m_opEBStencil
Definition: EBPoissonOp.H:292
static int getIndex(const IntVect &a_color)
Real m_alpha
Definition: EBPoissonOp.H:282
ProblemDomain m_domainCoarMG
Definition: EBPoissonOp.H:318
Piecewise constant interpolation.
Definition: EBMGInterp.H:33
const int SpaceDim
Definition: SPACE.H:38
virtual Real dotProduct(const LevelData< EBCellFAB > &a_1, const LevelData< EBCellFAB > &a_2)
VoF-centered stencil.
Definition: Stencils.H:60
EBMGInterp m_ebInterpMG
Definition: EBPoissonOp.H:315
virtual void applyOp(LevelData< EBCellFAB > &a_opPhi, const LevelData< EBCellFAB > &a_phi, bool a_homogeneousPhysBC, DataIterator &a_dit, bool do_exchange=true)
virtual void restrictResidual(LevelData< EBCellFAB > &a_resCoarse, LevelData< EBCellFAB > &a_phiFine, const LevelData< EBCellFAB > &a_rhsFine)
virtual void create(LevelData< EBCellFAB > &a_lhs, const LevelData< EBCellFAB > &a_rhs)
void levelJacobi(LevelData< EBCellFAB > &a_phi, const LevelData< EBCellFAB > &a_rhs)
RealVect m_dx
Definition: EBPoissonOp.H:277
virtual void setVal(LevelData< EBCellFAB > &a_lhs, const Real &a_value)
double Real
Definition: REAL.H:33
Definition: MultiGrid.H:30
LayoutData< BaseIVFAB< Real > > m_alphaDiagWeight
Definition: EBPoissonOp.H:298
void THIS_IS_AN_ERROR_MESSAGE(void)
Definition: EBPoissonOp.H:54
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:30
LayoutData< VoFIterator > m_vofItIrregDomLo[SpaceDim]
Definition: EBPoissonOp.H:306
LoHiSide
Definition: LoHiSide.H:27
virtual Real norm(const LevelData< EBCellFAB > &a_rhs, int a_ord)
int m_orderEB
Definition: EBPoissonOp.H:286
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
virtual void createCoarser(LevelData< EBCellFAB > &a_coarse, const LevelData< EBCellFAB > &a_fine, bool a_ghosted)
void applyOpRegularAllDirs(Box *a_loBox, Box *a_hiBox, int *a_hasLo, int *a_hasHi, Box &a_curOpPhiBox, Box &a_curPhiBox, int a_nComps, BaseFab< Real > &a_curOpPhiFAB, const BaseFab< Real > &a_curPhiFAB, bool a_homogeneousPhysBC, const DataIndex &a_dit, const Real &a_beta)
virtual void GSColorAllIrregular(LevelData< EBCellFAB > &a_phi, const LevelData< EBCellFAB > &a_rhs, const IntVect &a_color, const bool &a_homogeneousPhysBC, int icolor)
static IntVect getColor(const int &a_icolor)
RefCountedPtr< BaseDomainBC > m_domainBC
Definition: EBPoissonOp.H:274
LayoutData< BaseIVFAB< Real > > m_betaDiagWeight
Definition: EBPoissonOp.H:300
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
A Real vector in SpaceDim-dimensional space.
Definition: RealVect.H:41
void getInvDiagRHS(LevelData< EBCellFAB > &a_lhs, const LevelData< EBCellFAB > &a_rhs)
RefCountedPtr< BaseEBBC > m_ebBC
Definition: EBPoissonOp.H:275
Definition: DataIndex.H:114
Definition: EBPoissonOp.H:63
EBLevelGrid m_eblgCoarMG
Definition: EBPoissonOp.H:272
Piecewise constant interpolation.
Definition: EBMGAverage.H:31
LayoutData< RefCountedPtr< EBStencil > > m_rhsColorEBStencil[EBPO_NUMSTEN]
Definition: EBPoissonOp.H:296
LayoutData< VoFIterator > m_vofItIrregColorDomHi[EBPO_NUMSTEN][SpaceDim]
Definition: EBPoissonOp.H:310
void getOpVoFStencil(VoFStencil &a_stencil, const EBISBox &a_ebisbox, const VolIndex &a_vof)
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
bool m_hasMGObjects
Definition: EBPoissonOp.H:312
void getJacobiRelaxCoeff(LevelData< EBCellFAB > &a_relaxCoeff)
EBPoissonOp(const EBPoissonOp &a_opin)
Definition: EBPoissonOp.H:409
virtual void residual(LevelData< EBCellFAB > &a_residual, const LevelData< EBCellFAB > &a_phi, const LevelData< EBCellFAB > &a_rhs, bool a_homogeneousPhysBC=false)
Volume of Fluid Index.
Definition: VolIndex.H:31
void applyOpRegular(int idir, Box *a_loBox, Box *a_hiBox, int *a_hasLo, int *a_hasHi, Box &a_curOpPhiBox, Box &a_curPhiBox, int a_nComps, BaseFab< Real > &a_curOpPhiFAB, const BaseFab< Real > &a_curPhiFAB, bool a_homogeneousPhysBC, const DataIndex &a_dit, const Real &a_beta)
Definition: EBISLayout.H:39
virtual void relax(LevelData< EBCellFAB > &a_e, const LevelData< EBCellFAB > &a_residual, int a_iterations)
RealVect m_invDx
Definition: EBPoissonOp.H:279
virtual void prolongIncrement(LevelData< EBCellFAB > &a_phiThisLevel, const LevelData< EBCellFAB > &a_correctCoarse)
RealVect m_invDx2
Definition: EBPoissonOp.H:280
int m_relaxType
Definition: EBPoissonOp.H:288
Vector< IntVect > m_colors
Definition: EBPoissonOp.H:268
DisjointBoxLayout m_dblCoarMG
Definition: EBPoissonOp.H:316
EBMGAverage m_ebAverageMG
Definition: EBPoissonOp.H:314
EBLevelGrid m_eblg
Definition: EBPoissonOp.H:271
Real m_time
Definition: EBPoissonOp.H:284
virtual void assign(LevelData< EBCellFAB > &a_lhs, const LevelData< EBCellFAB > &a_rhs)