Chombo + EB  3.0
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 
291 
292  //stencils for operator evaluation
294  //stencils for operator evaluation on gauss-seidel colors
296  //stencils for lambda*rhs addition to opphi
298  //weights that get multiplied by alpha
300  //weights that get multiplied by beta
302 
303  //cache the irreg vofiterator
306 
309 
312 
314  //stuff below is only defined if m_hasMGObjects==true
320 
321  /// Get the coefficients (for all cells) to multiply the residual by
322  /// before adding it to the previous iterate for Jacobi iteration.
323  void getJacobiRelaxCoeff(LevelData<EBCellFAB>& a_relaxCoeff);
324 
325 private:
326 
327 
328  //internally useful but not for general consumption
329  //lots of hidden assumptions and all that
330 
331  void getOpVoFStencil(VoFStencil& a_stencil,
332  const EBISBox& a_ebisbox,
333  const VolIndex& a_vof);
334 
335  void getOpVoFStencil(VoFStencil& a_stencil,
336  const int& a_dir,
337  const Vector<VolIndex>& a_allMonotoneVoFs,
338  const EBISBox& a_ebisbox,
339  const VolIndex& a_vof,
340  const bool& a_lowOrder);
341 
342 
343  void getOpFaceStencil(VoFStencil& a_stencil,
344  const Vector<VolIndex>& a_allMonotoneVofs,
345  const EBISBox& a_ebisbox,
346  const VolIndex& a_vof,
347  int a_dir,
348  const Side::LoHiSide& a_side,
349  const FaceIndex& a_face,
350  const bool& a_lowOrder);
351 
352 
353  void levelJacobi(LevelData<EBCellFAB>& a_phi,
354  const LevelData<EBCellFAB>& a_rhs);
355 
356  ///
357  /** applyOp() on the regular cells. Only used in line relaxation stuff
358  opPhi comes in holding alpha*phi. this adds on beta*lapl phi*/
359  void applyOpRegular( int idir,
360  Box * a_loBox,
361  Box * a_hiBox,
362  int * a_hasLo,
363  int * a_hasHi,
364  Box & a_curOpPhiBox,
365  Box & a_curPhiBox,
366  int a_nComps,
367  BaseFab<Real> & a_curOpPhiFAB,
368  const BaseFab<Real> & a_curPhiFAB,
369  bool a_homogeneousPhysBC,
370  const DataIndex& a_dit,
371  const Real& a_beta);
372 
373  ///
374  /** applyOp() on the regular cells for all directions
375  opPhi comes in holding alpha*phi. this adds on beta*lapl phi*/
376  void applyOpRegularAllDirs(Box * a_loBox,
377  Box * a_hiBox,
378  int * a_hasLo,
379  int * a_hasHi,
380  Box & a_curOpPhiBox,
381  Box & a_curPhiBox,
382  int a_nComps,
383  BaseFab<Real> & a_curOpPhiFAB,
384  const BaseFab<Real> & a_curPhiFAB,
385  bool a_homogeneousPhysBC,
386  const DataIndex& a_dit,
387  const Real& a_beta);
388 
389  void applyDomainFlux(Box * a_loBox,
390  Box * a_hiBox,
391  int * a_hasLo,
392  int * a_hasHi,
393  Box & a_curPhiBox,
394  int a_nComps,
395  BaseFab<Real> & a_phiFAB,
396  bool a_homogeneousPhysBC,
397  const DataIndex& a_dit,
398  const Real& a_beta);
399 
400  ///
401  /**
402  This function computes: a_lhs = (1/diagonal)*a_rhs
403  Use this function to initialize the preconditioner
404  */
406  const LevelData<EBCellFAB>& a_rhs);
407 
408 private:
409  //copy constructor and operator= disallowed for all the usual reasons
410  EBPoissonOp(const EBPoissonOp& a_opin)
411  {
412  MayDay::Error("invalid operator");
413  }
414 
415  void operator=(const EBPoissonOp& a_opin)
416  {
417  MayDay::Error("invalid operator");
418  }
419 };
420 
421 
422 #include "NamespaceFooter.H"
423 #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:318
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:310
int m_numPreCondIters
Definition: EBPoissonOp.H:287
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:130
RealVect m_origin
Definition: EBPoissonOp.H:285
void operator=(const EBPoissonOp &a_opin)
Definition: EBPoissonOp.H:415
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: LayoutData.H:46
A strange but true thing to make copying from one boxlayoutdata to another fast.
Definition: Copier.H:137
LayoutData< RefCountedPtr< EBStencil > > m_colorEBStencil[EBPO_NUMSTEN]
Definition: EBPoissonOp.H:295
Definition: EBISBox.H:46
Definition: EBLevelGrid.H:30
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:304
Definition: DataIterator.H:140
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:308
LayoutData< VoFIterator > m_vofItIrregColor[EBPO_NUMSTEN]
Definition: EBPoissonOp.H:305
LayoutData< RefCountedPtr< EBStencil > > m_opEBStencil
Definition: EBPoissonOp.H:293
static int getIndex(const IntVect &a_color)
Real m_alpha
Definition: EBPoissonOp.H:282
ProblemDomain m_domainCoarMG
Definition: EBPoissonOp.H:319
Piecewise constant interpolation.
Definition: EBMGInterp.H:33
const int SpaceDim
Definition: SPACE.H:39
virtual Real dotProduct(const LevelData< EBCellFAB > &a_1, const LevelData< EBCellFAB > &a_2)
VoF-centered stencil.
Definition: Stencils.H:59
EBMGInterp m_ebInterpMG
Definition: EBPoissonOp.H:316
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:299
void THIS_IS_AN_ERROR_MESSAGE(void)
Definition: EBPoissonOp.H:54
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:31
LayoutData< VoFIterator > m_vofItIrregDomLo[SpaceDim]
Definition: EBPoissonOp.H:307
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:301
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
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:112
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:297
LayoutData< VoFIterator > m_vofItIrregColorDomHi[EBPO_NUMSTEN][SpaceDim]
Definition: EBPoissonOp.H:311
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:313
DataIterator dataIterator() const
Definition: LayoutDataI.H:79
void getJacobiRelaxCoeff(LevelData< EBCellFAB > &a_relaxCoeff)
EBPoissonOp(const EBPoissonOp &a_opin)
Definition: EBPoissonOp.H:410
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
Copier m_exchangeCopier
Definition: EBPoissonOp.H:290
int m_relaxType
Definition: EBPoissonOp.H:288
Vector< IntVect > m_colors
Definition: EBPoissonOp.H:268
DisjointBoxLayout m_dblCoarMG
Definition: EBPoissonOp.H:317
EBMGAverage m_ebAverageMG
Definition: EBPoissonOp.H:315
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)