ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ConvertForProjection.cpp File Reference
#include <AMReX.H>
#include <AMReX_MultiFab.H>
#include <ERF_Utils.H>
Include dependency graph for ERF_ConvertForProjection.cpp:

Functions

void ConvertForProjection (const MultiFab &den_div, const MultiFab &den_mlt, MultiFab &xmom, MultiFab &ymom, MultiFab &zmom, const Box &domain, const Vector< BCRec > &domain_bcs_type_h)
 
void compute_influx_outflux (Array< MultiFab *, AMREX_SPACEDIM > &vels_vec, Array< MultiFab *, AMREX_SPACEDIM > &area_vec, const Geometry &geom, Real &influx, Real &outflux)
 
void correct_outflow (const Geometry &geom_lev, Array< MultiFab *, AMREX_SPACEDIM > &vels_vec, const Box &domain, const Real alpha_fcf)
 
void enforceInOutSolvability (int, Array< MultiFab *, AMREX_SPACEDIM > &vels_vec, Array< MultiFab *, AMREX_SPACEDIM > &area_vec, const Geometry &geom)
 

Function Documentation

◆ compute_influx_outflux()

void compute_influx_outflux ( Array< MultiFab *, AMREX_SPACEDIM > &  vels_vec,
Array< MultiFab *, AMREX_SPACEDIM > &  area_vec,
const Geometry &  geom,
Real &  influx,
Real &  outflux 
)
183 {
184  influx = zero, outflux = zero;
185 
186  const Box domain = geom.Domain();
187  const auto& domlo = lbound(domain);
188  const auto& domhi = ubound(domain);
189 
190  // Normal face area (of undistorted mesh)
191  const Real* a_dx = geom.CellSize();
192  const Real ds_x = a_dx[1];
193  const Real ds_y = a_dx[0];
194 
195  IntVect ngrow = {0,0,0};
196 
197  // X-dir
198  auto const& vel_x = vels_vec[0]->const_arrays();
199  auto const& area_x = area_vec[0]->const_arrays();
200  influx += ds_x *
201  ParReduce(TypeList<ReduceOpSum>{},
202  TypeList<Real>{},
203  *vels_vec[0], ngrow,
204  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
205  noexcept -> GpuTuple<Real>
206  {
207  if ( (i == domlo.x && vel_x[box_no](i,j,k) > zero) ||
208  (i == domhi.x+1 && vel_x[box_no](i,j,k) < zero) ) {
209  return { std::abs(vel_x[box_no](i,j,k)) * area_x[box_no](i,j,k) };
210  } else {
211  return { zero };
212  }
213  });
214 
215  outflux += ds_x *
216  ParReduce(TypeList<ReduceOpSum>{},
217  TypeList<Real>{},
218  *vels_vec[0], ngrow,
219  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
220  noexcept -> GpuTuple<Real>
221  {
222  if ( (i == domlo.x && vel_x[box_no](i,j,k) < zero) ||
223  (i == domhi.x+1 && vel_x[box_no](i,j,k) > zero) ) {
224  return { std::abs(vel_x[box_no](i,j,k)) * area_x[box_no](i,j,k) };
225  } else {
226  return { zero };
227  }
228  });
229 
230  // Y-dir
231  auto const& vel_y= vels_vec[1]->const_arrays();
232  auto const& area_y = area_vec[1]->const_arrays();
233  influx += ds_y *
234  ParReduce(TypeList<ReduceOpSum>{},
235  TypeList<Real>{},
236  *vels_vec[1], ngrow,
237  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
238  noexcept -> GpuTuple<Real>
239  {
240  if ( (j == domlo.y && vel_y[box_no](i,j,k) > zero) ||
241  (j == domhi.y+1 && vel_y[box_no](i,j,k) < zero) ) {
242  return { std::abs(vel_y[box_no](i,j,k)) * area_y[box_no](i,j,k) };
243  } else {
244  return { zero };
245  }
246  });
247 
248  outflux += ds_y *
249  ParReduce(TypeList<ReduceOpSum>{},
250  TypeList<Real>{},
251  *vels_vec[1], ngrow,
252  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
253  noexcept -> GpuTuple<Real>
254  {
255  if ( (j == domlo.y && vel_y[box_no](i,j,k) < zero) ||
256  (j == domhi.y+1 && vel_y[box_no](i,j,k) > zero) ) {
257  return { std::abs(vel_y[box_no](i,j,k)) * area_y[box_no](i,j,k) };
258  } else {
259  return { zero };
260  }
261  });
262 
263  ParallelDescriptor::ReduceRealSum(influx);
264  ParallelDescriptor::ReduceRealSum(outflux);
265 }
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19

Referenced by enforceInOutSolvability().

Here is the caller graph for this function:

◆ ConvertForProjection()

void ConvertForProjection ( const MultiFab &  den_div,
const MultiFab &  den_mlt,
MultiFab &  xmom,
MultiFab &  ymom,
MultiFab &  zmom,
const Box &  domain,
const Vector< BCRec > &  domain_bcs_type_h 
)

Convert momentum to velocity by dividing by density averaged onto faces

Parameters
[out]xvelx-component of velocity
[out]yvely-component of velocity
[out]zvelz-component of velocity
[in]densitydensity at cell centers
[in]xmom_inx-component of momentum
[in]ymom_iny-component of momentum
[in]zmom_inz-component of momentum
[in]domainDomain at this level
[in]domain_bcs_type_hhost vector for domain boundary conditions
28 {
29  BL_PROFILE_VAR("ConvertForProjection()",onvertForProjection);
30 
31  const BCRec* bc_ptr_h = domain_bcs_type_h.data();
32 
33 #ifdef _OPENMP
34 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
35 #endif
36  for ( MFIter mfi(den_div,TilingIfNotGPU()); mfi.isValid(); ++mfi)
37  {
38  // We need velocity in the interior ghost cells (init == real)
39  Box bx = mfi.tilebox();
40 
41  Box tbx = surroundingNodes(bx,0);
42  Box tby = surroundingNodes(bx,1);
43  Box tbz = surroundingNodes(bx,2);
44 
45  if ( (bx.smallEnd(0) == domain.smallEnd(0)) &&
46  ( (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir) ||
47  (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir_upwind) ) )
48  {
49  tbx.growLo(0,-1);
50  }
51  if ( (bx.bigEnd(0) == domain.bigEnd(0)) &&
52  ( (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir) ||
53  (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir_upwind) ) )
54  {
55  tbx.growHi(0,-1);
56  }
57  if ( (bx.smallEnd(1) == domain.smallEnd(1)) &&
58  ( (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir) ||
59  (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir_upwind) ) )
60  {
61  tby.growLo(1,-1);
62  }
63  if ( (bx.bigEnd(1) == domain.bigEnd(1)) &&
64  ( (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir) ||
65  (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir_upwind) ) )
66  {
67  tby.growHi(1,-1);
68  }
69 
70  // Conserved variables on cell centers -- we use this for density
71  const Array4<const Real>& den_div_arr = den_div.const_array(mfi);
72  const Array4<const Real>& den_mlt_arr = den_mlt.const_array(mfi);
73 
74  // Momentum on faces
75  Array4<Real> const& momx = xmom.array(mfi);
76  Array4<Real> const& momy = ymom.array(mfi);
77  Array4<Real> const& momz = zmom.array(mfi);
78 
79  ParallelFor(tbx, tby, tbz,
80  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
81  momx(i,j,k) *= ( den_mlt_arr(i,j,k,Rho_comp) + den_mlt_arr(i-1,j,k,Rho_comp) )
82  / ( den_div_arr(i,j,k,Rho_comp) + den_div_arr(i-1,j,k,Rho_comp) );
83 
84  },
85  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
86  momy(i,j,k) *= ( den_mlt_arr(i,j,k,Rho_comp) + den_mlt_arr(i,j-1,k,Rho_comp) )
87  / ( den_div_arr(i,j,k,Rho_comp) + den_div_arr(i,j-1,k,Rho_comp) );
88  },
89  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
90  momz(i,j,k) *= ( den_mlt_arr(i,j,k,Rho_comp) + den_mlt_arr(i,j,k-1,Rho_comp) )
91  / ( den_div_arr(i,j,k,Rho_comp) + den_div_arr(i,j,k-1,Rho_comp) );
92  });
93 
94  if (bx.smallEnd(0) == domain.smallEnd(0)) {
95  if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir)
96  {
97  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
98  momx(i,j,k) *= den_mlt_arr(i-1,j,k,Rho_comp) / den_div_arr(i-1,j,k,Rho_comp) ;
99  });
100  }
101  else if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir_upwind)
102  {
103  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
104  if (momx(i,j,k) >= zero) {
105  momx(i,j,k) *= den_mlt_arr(i-1,j,k,Rho_comp) / den_div_arr(i-1,j,k,Rho_comp) ;
106  } else {
107  momx(i,j,k) *= ( den_mlt_arr(i,j,k,Rho_comp) + den_mlt_arr(i-1,j,k,Rho_comp) )
108  / ( den_div_arr(i,j,k,Rho_comp) + den_div_arr(i-1,j,k,Rho_comp) );
109  }
110  });
111  }
112  }
113 
114  if (bx.bigEnd(0) == domain.bigEnd(0)) {
115  if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir)
116  {
117  ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
118  momx(i,j,k) *= den_mlt_arr(i,j,k,Rho_comp) / den_div_arr(i,j,k,Rho_comp) ;
119  });
120  }
121  else if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir_upwind)
122  {
123  ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
124  if (momx(i,j,k) <= zero) {
125  momx(i,j,k) *= den_mlt_arr(i,j,k,Rho_comp) / den_div_arr(i,j,k,Rho_comp) ;
126  } else {
127  momx(i,j,k) *= ( den_mlt_arr(i,j,k,Rho_comp) + den_mlt_arr(i-1,j,k,Rho_comp) )
128  / ( den_div_arr(i,j,k,Rho_comp) + den_div_arr(i-1,j,k,Rho_comp) );
129  }
130  });
131  }
132  }
133 
134  if (bx.smallEnd(1) == domain.smallEnd(1)) {
135  if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir)
136  {
137  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
138  momy(i,j,k) *= den_mlt_arr(i,j-1,k,Rho_comp) / den_div_arr(i,j-1,k,Rho_comp) ;
139  });
140  }
141  else if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir_upwind)
142  {
143  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
144  if (momy(i,j,k) >= zero) {
145  momy(i,j,k) *= den_mlt_arr(i,j-1,k,Rho_comp) / den_div_arr(i,j-1,k,Rho_comp) ;
146  } else {
147  momy(i,j,k) *= ( den_mlt_arr(i,j,k,Rho_comp) + den_mlt_arr(i,j-1,k,Rho_comp) )
148  / ( den_div_arr(i,j,k,Rho_comp) + den_div_arr(i,j-1,k,Rho_comp) );
149  }
150  });
151  }
152  }
153 
154  if (bx.bigEnd(1) == domain.bigEnd(1)) {
155  if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir)
156  {
157  ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
158  momy(i,j,k) *= den_mlt_arr(i,j,k,Rho_comp) / den_div_arr(i,j,k,Rho_comp) ;
159  });
160  }
161  else if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir_upwind)
162  {
163  ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
164  if (momy(i,j,k) <= zero) {
165  momy(i,j,k) *= den_mlt_arr(i,j,k,Rho_comp) / den_div_arr(i,j,k,Rho_comp) ;
166  } else {
167  momy(i,j,k) *= ( den_mlt_arr(i,j,k,Rho_comp) + den_mlt_arr(i,j-1,k,Rho_comp) )
168  / ( den_div_arr(i,j,k,Rho_comp) + den_div_arr(i,j-1,k,Rho_comp) );
169  }
170  });
171  }
172  }
173 
174 
175  } // end MFIter
176 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const Real *dx=geomdata.CellSize();const Real x=(i+0.5) *dx[0];const Real y=(j+0.5) *dx[1];const Real Omg=erf_vortex_Gaussian(x, y, xc, yc, R, beta, sigma);const Real deltaT=-(gamma - 1.0)/(2.0 *sigma *sigma) *Omg *Omg;const Real rho_norm=std::pow(1.0+deltaT, inv_gm1);const Real T=(1.0+deltaT) *T_inf;const Real p=std::pow(rho_norm, Gamma)/Gamma *rho_0 *a_inf *a_inf;const Real rho_theta=rho_0 *rho_norm *(T *std::pow(p_0/p, rdOcp));state_pert(i, j, k, RhoTheta_comp)=rho_theta - getRhoThetagivenP(p_hse(i, j, k));const Real r2d_xy=std::sqrt((x-xc) *(x-xc)+(y-yc) *(y-yc));state_pert(i, j, k, RhoScalar_comp)=0.25 *(1.0+std::cos(PI *std::min(r2d_xy, R)/R));})
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:76
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:227
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:235
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:178
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:179
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:177

Referenced by ERF::project_momenta().

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

◆ correct_outflow()

void correct_outflow ( const Geometry &  geom_lev,
Array< MultiFab *, AMREX_SPACEDIM > &  vels_vec,
const Box &  domain,
const Real  alpha_fcf 
)
272 {
273  for (OrientationIter oit; oit != nullptr; ++oit) {
274  const auto ori = oit();
275  const int dir = ori.coordDir();
276  const auto oriIsLow = ori.isLow();
277  const auto oriIsHigh = ori.isHigh();
278 
279  if (dir < 2) {
280 
281  // MultiFab for normal velocity
282  const auto& vel_mf = vels_vec[dir];
283 
284  // Domain extent indices for the velocities
285  int dlo = domain.smallEnd(dir);
286  int dhi = domain.bigEnd(dir)+1;
287 
288  // get BCs for the normal velocity and set the boundary index
289  int bndry;
290  if (oriIsLow) {
291  bndry = dlo;
292  } else {
293  bndry = dhi;
294  }
295 
296  // Assume here that all domain boundaries are viewed as direction_dependent!
297  if (!geom_lev.isPeriodic(dir))
298  {
299  for (MFIter mfi(*vel_mf, false); mfi.isValid(); ++mfi) {
300 
301  Box box = mfi.validbox();
302 
303  // Enter further only if the box boundary is at the domain boundary
304  if ((oriIsLow && (box.smallEnd(dir) == dlo))
305  || (oriIsHigh && (box.bigEnd(dir) == dhi))) {
306 
307  // create a 2D box normal to dir at the low/high boundary
308  Box box2d(box); box2d.setRange(dir, bndry);
309 
310  auto vel_arr = vel_mf->array(mfi);
311 
312  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
313  {
314  if ((oriIsLow && vel_arr(i,j,k) < 0)
315  || (oriIsHigh && vel_arr(i,j,k) > 0)) {
316  vel_arr(i,j,k) *= alpha_fcf;
317  }
318  });
319  }
320  } // mfi
321  } // geom
322  } // dir < 2
323  } // ori
324 }

Referenced by enforceInOutSolvability().

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

◆ enforceInOutSolvability()

void enforceInOutSolvability ( int  ,
Array< MultiFab *, AMREX_SPACEDIM > &  vels_vec,
Array< MultiFab *, AMREX_SPACEDIM > &  area_vec,
const Geometry &  geom 
)
330 {
331  Real small_vel = Real(1.e-8);
332 
333  Real influx = zero, outflux = zero;
334 
335  const Box domain = geom.Domain();
336 
337  Real influx_lev = zero, outflux_lev = zero;
338  compute_influx_outflux(vels_vec, area_vec, geom, influx_lev, outflux_lev);
339  influx += influx_lev;
340  outflux += outflux_lev;
341  amrex::Print() <<" TOTAL INFLUX / OUTFLOW " << influx << " " << outflux << std::endl;
342 
343  if ((influx > small_vel) && (outflux < small_vel)) {
344  Abort("Cannot enforce solvability, no outflow from the direction dependent boundaries");
345  } else if ((influx < small_vel) && (outflux < small_vel)) {
346  return; // do nothing
347  } else {
348  const Real alpha_fcf = influx/outflux; // flux correction factor
349  correct_outflow(geom, vels_vec, domain, alpha_fcf);
350 
351  // Just for diagnostic purposes!
352  // Real influx_lev = zero, outflux_lev = zero;
353  // compute_influx_outflux(vels_vec, area_vec, geom, influx_lev, outflux_lev);
354  // amrex::Print() <<" TOTAL INFLUX / OUTFLOW " << influx_lev << " " << outflux_lev << std::endl;
355  }
356 }
correct_outflow
void correct_outflow(const Geometry &geom_lev, Array< MultiFab *, AMREX_SPACEDIM > &vels_vec, const Box &domain, const Real alpha_fcf)
Definition: ERF_ConvertForProjection.cpp:267
compute_influx_outflux
void compute_influx_outflux(Array< MultiFab *, AMREX_SPACEDIM > &vels_vec, Array< MultiFab *, AMREX_SPACEDIM > &area_vec, const Geometry &geom, Real &influx, Real &outflux)
Definition: ERF_ConvertForProjection.cpp:178

Referenced by ERF::project_momenta().

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