ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Private Attributes | List of all members
eb_aux_ Class Reference

#include <ERF_EBAux.H>

Collaboration diagram for eb_aux_:
Collaboration graph
[legend]

Public Member Functions

 eb_aux_ ()
 
 ~eb_aux_ ()
 
void define (int const &a_idim, amrex::Geometry const &a_geom, amrex::BoxArray const &a_grids, amrex::DistributionMapping const &a_dmap, amrex::Vector< int > const &a_ngrow, amrex::EBFArrayBoxFactory const *a_factory)
 
void set_verbose ()
 

Private Attributes

int m_verbose
 
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags = nullptr
 
amrex::MultiFab * m_volfrac = nullptr
 

Constructor & Destructor Documentation

◆ eb_aux_()

eb_aux_::eb_aux_ ( )
13  : m_verbose(0)
14 // ,m_defined(0)
15 {}
eb_aux_::m_verbose
int m_verbose
Definition: ERF_EBAux.H:29

◆ ~eb_aux_()

eb_aux_::~eb_aux_ ( )
8 {
9 }

Member Function Documentation

◆ define()

void eb_aux_::define ( int const &  a_idim,
amrex::Geometry const &  a_geom,
amrex::BoxArray const &  a_grids,
amrex::DistributionMapping const &  a_dmap,
amrex::Vector< int > const &  a_ngrow,
amrex::EBFArrayBoxFactory const *  a_factory 
)
25 {
26  // Box dbox(a_geom.Domain());
27 
28  const IntVect vdim(IntVect::TheDimensionVector(a_idim));
29 
30  const BoxArray& grids = amrex::convert(a_grids, vdim);
31 
32  m_cellflags = new FabArray<EBCellFlagFab>(grids, a_dmap, 1, a_ngrow[0], MFInfo(),
33  DefaultFabFactory<EBCellFlagFab>());
34 
35  m_cellflags->setVal(EBCellFlag::TheDefaultCell());
36 
37  m_volfrac = new MultiFab(grids, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory());
38 
39 
40 #if 0
41  m_centroid = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, m_ngrow[1], *m_cellflags);
42  m_bndrycent = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, m_grow[2], *m_cellflags);
43 
44  m_bndryarea = new MultiCutFab(a_ba, a_dm, 1, m_grow[2], *m_cellflags);
45  m_bndrynorm = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, m_grow[2], *m_cellflags);
46 
47  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
48  const BoxArray& faceba = amrex::convert(a_ba, IntVect::TheDimensionVector(idim));
49  m_areafrac[idim] = new MultiCutFab(faceba, a_dm, 1, m_grow[2], *m_cellflags);
50  m_facecent[idim] = new MultiCutFab(faceba, a_dm, AMREX_SPACEDIM-1, m_grow[2], *m_cellflags);
51  }
52 #endif
53 
54  const auto& FlagFab = a_factory->getMultiEBCellFlagFab();
55 
56  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
57 
58  const Box& bx = mfi.validbox();
59 
60  if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
61 
62  Array4<EBCellFlag const> const& flag = FlagFab.const_array(mfi);
63 
64  Array4<Real const> const& vfrac = (a_factory->getVolFrac()).const_array(mfi);
65  // Array4<Real const> const& ccent = (a_factory->getCentroid()).const_array(mfi);
66 
67  Array4<Real const> const& afrac = (a_factory->getAreaFrac()[a_idim])->const_array(mfi);
68 
69  // EB normal and face centroid
70  Array4<Real const> const& bnorm = a_factory->getBndryNormal()[mfi].const_array();
71  Array4<Real const> const& bcent = a_factory->getBndryCent()[mfi].const_array();
72 
73  Array4<EBCellFlag> const& aux_flag = m_cellflags->array(mfi);
74  Array4<Real> const& aux_vfrac = m_volfrac->array(mfi);
75 
76  ParallelFor(bx, [
77 #ifndef AMREX_USE_GPU
78  verbose=m_verbose,
79 #endif
80  vfrac, afrac, bnorm, bcent, flag,
81  aux_flag, aux_vfrac, vdim, idim=a_idim ]
82  AMREX_GPU_DEVICE (int i, int j, int k) noexcept
83  {
84  aux_flag(i,j,k).setCovered();
85  aux_flag(i,j,k).setDisconnected();
86 
87  aux_vfrac(i,j,k) = 0.0;
88 
89  IntVect const iv_hi(i,j,k);
90  IntVect const iv_lo(iv_hi - vdim);
91 
92  if ( flag(iv_lo).isRegular() && flag(iv_hi).isRegular()) {
93 
94  aux_flag(i,j,k).setRegular();
95  aux_flag(i,j,k).setConnected(vdim);
96 
97  aux_vfrac(i,j,k) = 1.0;
98 
99  } else if ( flag(iv_lo).isCovered() && flag(iv_hi).isCovered()) {
100 
101  // defaults to covered and disconnected.
102 
103  } else {
104 
105 
106 #ifndef AMREX_USE_GPU
107  if (verbose) { Print() << "\ncell: " << amrex::IntVect(i,j,k) << "\n"; }
108 #endif
109  Array<Real,AMREX_SPACEDIM> lo_arr = {-0.5,-0.5,-0.5};
110  Array<Real,AMREX_SPACEDIM> hi_arr = { 0.5, 0.5, 0.5};
111 
112  // plane point and normal
113  RealVect lo_point(bcent(iv_lo,0), bcent(iv_lo,1), bcent(iv_lo,2));
114  RealVect lo_normal(bnorm(iv_lo,0), bnorm(iv_lo,1), bnorm(iv_lo,2));
115 
116  // High side of low cell
117  lo_arr[idim] = 0.0;
118  hi_arr[idim] = 0.5;
119  RealBox lo_rbx(lo_arr.data(), hi_arr.data());
120 
121  eb_cut_cell_ lo_eb_cc(flag(iv_lo), lo_rbx, lo_point, lo_normal);
122 
123  // cell iv_lo covered (regular) imples lo_eb_cc is covered (regular)
124  // The inverse is not always true.
125  AMREX_ASSERT( !flag(iv_lo).isCovered() || lo_eb_cc.isCovered() );
126  AMREX_ASSERT( !flag(iv_lo).isRegular() || lo_eb_cc.isRegular() );
127 
128  // plane point and normal
129  RealVect hi_point(bcent(iv_hi,0), bcent(iv_hi,1), bcent(iv_hi,2));
130  RealVect hi_normal(bnorm(iv_hi,0), bnorm(iv_hi,1), bnorm(iv_hi,2));
131 
132  // Low side of high cell
133  lo_arr[idim] = -0.5;
134  hi_arr[idim] = 0.0;
135  RealBox hi_rbx(lo_arr.data(), hi_arr.data());
136 
137  eb_cut_cell_ hi_eb_cc(flag(iv_hi), hi_rbx, hi_point, hi_normal);
138 
139  // cell iv_hi covered (regular) imples hi_eb_cc is covered (regular)
140  // The inverse is not always true.
141  AMREX_ASSERT( !flag(iv_hi).isCovered() || hi_eb_cc.isCovered() );
142  AMREX_ASSERT( !flag(iv_hi).isRegular() || hi_eb_cc.isRegular() );
143 
144 #if defined(AMREX_DEBUG) || defined(AMREX_TESTING) || 1
145 
146  { /***************************** SANITY CHECK ***********************\
147  * Perform some basic sanity checks to verify that what we computed *
148  * for cell (i,j,k) compares to what we know to be true. *
149  \******************************************************************/
150 
151  // Compute the cut-cell for the high side of the high cell. This is
152  // only needed for sanity checks.
153 
154  eb_cut_cell_ hi_hi_eb_cc(flag(iv_hi), lo_rbx, hi_point, hi_normal);
155 
156  // cell iv_hi covered (regular) imples hi_hi_eb_cc is covered (regular)
157  // The inverse is not always true.
158 #ifndef AMREX_USE_GPU
159  if ( !(!flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular()) ||
160  !(!flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered()) ) {
161  Print() << "flag(iv_hi) and hi_hi_eb_cc flags do not agree\n"
162  << "\n isRegular() " << flag(iv_hi).isRegular() << " " << hi_hi_eb_cc.isRegular()
163  << "\n isCovered() " << flag(iv_hi).isCovered() << " " << hi_hi_eb_cc.isCovered()
164  << "\n";
165  }
166 #endif
167  // If cell iv_hi is regular or covered, then hi_hi_eb_cc must also
168  // be regular or covered. The inverse is not true.
169  AMREX_ALWAYS_ASSERT( !flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular() );
170  AMREX_ALWAYS_ASSERT( !flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered() );
171 
172  // The area and volume fractions that are computed for the scalar grid
173  // are slightly different than those we compute from the geometric
174  // reconstruction using the EB point and normal. However, we expect
175  // that the area fractions computed here will give back the same
176  // normal we used to compute them.
177  if ( flag(iv_hi).isSingleValued() ) {
178 
179  Real const adx = (idim == 0)
180  ? (hi_eb_cc.areaLo(0) - hi_hi_eb_cc.areaHi(0))
181  : (hi_eb_cc.areaLo(0) + hi_hi_eb_cc.areaLo(0))
182  - (hi_eb_cc.areaHi(0) + hi_hi_eb_cc.areaHi(0));
183 
184  Real const ady = (idim == 1)
185  ? (hi_eb_cc.areaLo(1) - hi_hi_eb_cc.areaHi(1))
186  : (hi_eb_cc.areaLo(1) + hi_hi_eb_cc.areaLo(1))
187  - (hi_eb_cc.areaHi(1) + hi_hi_eb_cc.areaHi(1));
188 
189  Real const adz = (idim == 2)
190  ? (hi_eb_cc.areaLo(2) - hi_hi_eb_cc.areaHi(2))
191  : (hi_eb_cc.areaLo(2) + hi_hi_eb_cc.areaLo(2))
192  - (hi_eb_cc.areaHi(2) + hi_hi_eb_cc.areaHi(2));
193 
194  Real const apnorm = std::sqrt(adx*adx + ady*ady + adz*adz);
195 
196  // EB normal
197  Real const apnorminv = 1. / apnorm;
198  RealVect const normal(adx*apnorminv, ady*apnorminv, adz*apnorminv);
199  Real const dot_normals = normal.dotProduct(hi_normal);
200 
201 #ifndef AMREX_USE_GPU
202  if ( !amrex::almostEqual(dot_normals, 1.0) ) {
203  Print() << "\nFail: check-1 dot_normals " << dot_normals
204  << '\n';
205 
206  hi_eb_cc.debug();
207  hi_hi_eb_cc.debug();
208 
209  } else if (verbose) {
210  Print() << "Pass: dot_normals = 1.0\n";
211 
212  }
213 #endif
214  AMREX_ALWAYS_ASSERT( amrex::almostEqual(dot_normals, 1.0) );
215  }
216 
217  // The idim area of hi_eb_cc.areaHi() should equal hi_hi_eb_cc.areaLo()
218  {
219 #ifndef AMREX_USE_GPU
220  Real const abs_err = std::abs( hi_eb_cc.areaHi(idim) - hi_hi_eb_cc.areaLo(idim) );
221  Real machine_tol = 10.0*std::numeric_limits<amrex::Real>::epsilon();
222  if ( abs_err >= machine_tol ) {
223  Print() << "\nFail: check-4 area abs_err: " << abs_err
224  << "\n hi_eb_cc.areaHi " << hi_eb_cc.areaHi(idim)
225  << "\n hi_hi_eb_cc.areaLo " << hi_hi_eb_cc.areaLo(idim)
226  << '\n';
227  } else if (verbose) {
228  Print() << "Pass: hi_eb_cc.areaHi = hi_hi_eb_cc.areaLo"
229  << " abs_err: " << abs_err << "\n";
230  }
231  AMREX_ALWAYS_ASSERT( abs_err < machine_tol );
232 #endif
233  }
234 
235  // The low-side area of hi_eb_cc should equal idim afrac.
236  { Real const abs_err = amrex::min(std::abs(lo_eb_cc.areaHi(idim) - afrac(iv_hi)),
237  std::abs(hi_eb_cc.areaLo(idim) - afrac(iv_hi)));
238  Real compare_tol = 5.0e-6;
239 #ifndef AMREX_USE_GPU
240  if ( abs_err >= compare_tol ) {
241  //hi_eb_cc.debug();
242  Print() << "\nFail: check-2 area abs_err " << abs_err
243  << "\n hi_eb_cc.areaLo(" << idim << ") = " << hi_eb_cc.areaLo(idim)
244  << "\n lo_eb_cc.areaHi(" << idim << ") = " << lo_eb_cc.areaHi(idim)
245  << "\n afrac" << iv_hi << " = " << afrac(iv_hi)
246  << '\n';
247  } else if (verbose) {
248  Print() << "Pass: hi_eb_cc.areaLo = afrac = " << afrac(iv_hi)
249  << " abs_err: " << abs_err << "\n";
250  }
251 #endif
252  AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
253  }
254 
255  // The combined volumes of hi_eb_cc.areaHi() and hi_hi_eb_cc should
256  // equal vfrac(iv_hi).
257  { Real const vol = hi_eb_cc.volume() + hi_hi_eb_cc.volume();
258  Real const abs_err = amrex::Math::abs(vfrac(iv_hi) - vol);
259  Real compare_tol = 5.0e-6;
260 #ifndef AMREX_USE_GPU
261  if ( abs_err >= compare_tol ) {
262  hi_eb_cc.debug();
263  hi_hi_eb_cc.debug();
264  amrex::Print() << "\nFail: check-4 volume abs_err: " << abs_err
265  << "\n point: " << hi_point
266  << "\n normal: " << hi_normal
267  << "\n hi_eb_cc.volume() " << hi_eb_cc.volume()
268  << "\n hi_hi_eb_cc.volume() " << hi_hi_eb_cc.volume()
269  << "\n vfrac: " << vfrac(iv_hi)
270  << '\n';
271  } else if (verbose) {
272  Print() << "Pass: hi_eb_cc + hi_hi_eb_cc = vfrac = " << vfrac(iv_hi)
273  << " abs_err: " << abs_err << "\n";
274  }
275 #endif
276  AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
277  }
278  }
279 #endif
280 
281  if (lo_eb_cc.isCovered() && hi_eb_cc.isCovered()) {
282 
283  // defaults to covered and disconnected.
284 
285  } else if (lo_eb_cc.isRegular() && hi_eb_cc.isRegular()) {
286 
287  aux_flag(i,j,k).setRegular();
288  aux_flag(i,j,k).setConnected(vdim);
289 
290  aux_vfrac(i,j,k) = 1.0;
291 
292  } else {
293 
294  if (lo_eb_cc.isCovered()) { }
295  if (hi_eb_cc.isCovered()) { }
296 
297  }
298  }
299  });
300 
301  }
302 
303  }
304 }
eb_aux_::m_cellflags
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags
Definition: ERF_EBAux.H:33
eb_aux_::m_volfrac
amrex::MultiFab * m_volfrac
Definition: ERF_EBAux.H:35
eb_cut_cell_
Definition: ERF_EBCutCell.H:12

Referenced by eb_::make_factory().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_verbose()

void eb_aux_::set_verbose ( )
inline
25 { m_verbose = 1; }

Member Data Documentation

◆ m_cellflags

amrex::FabArray<amrex::EBCellFlagFab>* eb_aux_::m_cellflags = nullptr
private

Referenced by define().

◆ m_verbose

int eb_aux_::m_verbose
private

Referenced by define(), and set_verbose().

◆ m_volfrac

amrex::MultiFab* eb_aux_::m_volfrac = nullptr
private

Referenced by define().

The documentation for this class was generated from the following files: