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

Functions

void compute_influx_outflux_bdy (FArrayBox &bdy_data_xlo, FArrayBox &bdy_data_xhi, FArrayBox &bdy_data_ylo, FArrayBox &bdy_data_yhi, Array< MultiFab *, AMREX_SPACEDIM > &area_vec, const Geometry &geom, Real &influx, Real &outflux, const int n)
 
void enforceInOutSolvability_bdy (const MultiFab &rho0, FArrayBox &bdy_data_xlo, FArrayBox &bdy_data_xhi, FArrayBox &bdy_data_ylo, FArrayBox &bdy_data_yhi, Array< MultiFab *, AMREX_SPACEDIM > &area_vec, const Geometry &geom, const Vector< BCRec > &domain_bcs_type_h)
 

Function Documentation

◆ compute_influx_outflux_bdy()

void compute_influx_outflux_bdy ( FArrayBox &  bdy_data_xlo,
FArrayBox &  bdy_data_xhi,
FArrayBox &  bdy_data_ylo,
FArrayBox &  bdy_data_yhi,
Array< MultiFab *, AMREX_SPACEDIM > &  area_vec,
const Geometry &  geom,
Real &  influx,
Real &  outflux,
const int  n 
)
158 {
159  BL_PROFILE_VAR("compute_influx_outflux_bdy()", computeInfluxOutfluxBdy);
160 
161  influx = zero, outflux = zero;
162 
163  const Box domain = geom.Domain();
164  const auto& domlo = lbound(domain);
165  const auto& domhi = ubound(domain);
166 
167  // Normal face area (of undistorted mesh)
168  const Real* a_dx = geom.CellSize();
169  const Real ds_x = a_dx[1];
170  const Real ds_y = a_dx[0];
171 
172  IntVect ngrow = {0,0,0};
173 
174  // X-dir
175  const auto& bdatxlo = bdy_data_xlo.const_array();
176  const auto& bdatxhi = bdy_data_xhi.const_array();
177  auto const& area_x = area_vec[0]->const_arrays();
178  influx += ds_x *
179  ParReduce(TypeList<ReduceOpSum>{},
180  TypeList<Real>{},
181  *area_vec[0], ngrow,
182  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
183  noexcept -> GpuTuple<Real>
184  {
185  if ( (i == domlo.x + n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxlo(i,j,k) > zero) {
186  return { std::abs(bdatxlo(i,j,k)) * area_x[box_no](i,j,k) };
187  } else if ( (i == domhi.x+1 - n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxhi(i,j,k) < zero) {
188  return { std::abs(bdatxhi(i,j,k)) * area_x[box_no](i,j,k) };
189  } else {
190  return { zero };
191  }
192  });
193 
194  outflux += ds_x *
195  ParReduce(TypeList<ReduceOpSum>{},
196  TypeList<Real>{},
197  *area_vec[0], ngrow,
198  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
199  noexcept -> GpuTuple<Real>
200  {
201  if (i == (domlo.x + n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxlo(i,j,k) < zero) {
202  return { std::abs(bdatxlo(i,j,k)) * area_x[box_no](i,j,k) };
203  } else if ( (i == domhi.x+1 - n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxhi(i,j,k) > zero) {
204  return { std::abs(bdatxhi(i,j,k)) * area_x[box_no](i,j,k) };
205  } else {
206  return { zero };
207  }
208  });
209 
210  // Y-dir
211  const auto& bdatylo = bdy_data_ylo.const_array();
212  const auto& bdatyhi = bdy_data_yhi.const_array();
213  auto const& area_y = area_vec[1]->const_arrays();
214  influx += ds_y *
215  ParReduce(TypeList<ReduceOpSum>{},
216  TypeList<Real>{},
217  *area_vec[1], ngrow,
218  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
219  noexcept -> GpuTuple<Real>
220  {
221  if ( (j == domlo.y + n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatylo(i,j,k) > zero) {
222  return { std::abs(bdatylo(i,j,k)) * area_y[box_no](i,j,k) };
223  } else if ( (j == domhi.y+1 - n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatyhi(i,j,k) < zero) {
224  return { std::abs(bdatyhi(i,j,k)) * area_y[box_no](i,j,k) };
225  } else {
226  return { zero };
227  }
228  });
229 
230  outflux += ds_y *
231  ParReduce(TypeList<ReduceOpSum>{},
232  TypeList<Real>{},
233  *area_vec[1], ngrow,
234  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
235  noexcept -> GpuTuple<Real>
236  {
237  if ( (j == domlo.y + n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatylo(i,j,k) < zero) {
238  return { std::abs(bdatylo(i,j,k)) * area_y[box_no](i,j,k) };
239  } else if ( (j == domhi.y+1 - n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatyhi(i,j,k) > zero) {
240  return { std::abs(bdatyhi(i,j,k)) * area_y[box_no](i,j,k) };
241  } else {
242  return { zero };
243  }
244  });
245 
246  ParallelDescriptor::ReduceRealSum(influx);
247  ParallelDescriptor::ReduceRealSum(outflux);
248 }
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19

Referenced by enforceInOutSolvability_bdy().

Here is the caller graph for this function:

◆ enforceInOutSolvability_bdy()

void enforceInOutSolvability_bdy ( const MultiFab &  rho0,
FArrayBox &  bdy_data_xlo,
FArrayBox &  bdy_data_xhi,
FArrayBox &  bdy_data_ylo,
FArrayBox &  bdy_data_yhi,
Array< MultiFab *, AMREX_SPACEDIM > &  area_vec,
const Geometry &  geom,
const Vector< BCRec > &  domain_bcs_type_h 
)
258 {
259  BL_PROFILE_VAR("enforceInOutSolvability_bdy()", enforceInOutSolvabilityBdy);
260 
261 #ifdef AMREX_USE_FLOAT
262  Real small_vel = Real(1.e-6);
263 #else
264  Real small_vel = Real(1.e-8);
265 #endif
266 
267  Box domain(geom.Domain());
268 
269  bool multiply_by_rho0 = true;
270  scale_bdy_normal_by_rho0(rho0, bdy_data_xlo, 0, domain, domain_bcs_type_h, multiply_by_rho0);
271  scale_bdy_normal_by_rho0(rho0, bdy_data_xhi, 0, domain, domain_bcs_type_h, multiply_by_rho0);
272  scale_bdy_normal_by_rho0(rho0, bdy_data_ylo, 1, domain, domain_bcs_type_h, multiply_by_rho0);
273  scale_bdy_normal_by_rho0(rho0, bdy_data_yhi, 1, domain, domain_bcs_type_h, multiply_by_rho0);
274 
275  int width =bdy_data_xlo.box().length(0);
276  AMREX_ASSERT(width == bdy_data_xhi.box().length(0));
277  AMREX_ASSERT(width == bdy_data_ylo.box().length(1));
278  AMREX_ASSERT(width == bdy_data_yhi.box().length(1));
279  // amrex::Print() << "width in enforce solvability " << width << std::endl;
280 
281  for (int n = 0; n < width; n++)
282  {
283  Real influx = zero, outflux = zero;
284  compute_influx_outflux_bdy(bdy_data_xlo, bdy_data_xhi,
285  bdy_data_ylo, bdy_data_yhi,
286  area_vec, geom, influx, outflux, n);
287 
288  // amrex::Print() << " TOTAL BDY INFLUX / OUTFLOW AT WIDTH " << n << " " << influx << " " << outflux << "\n";
289 
290  if ((influx > small_vel) && (outflux < small_vel)) {
291  Abort("Cannot enforce solvability on bdy_data, no outflow from the direction dependent boundaries");
292  } else if ((influx < small_vel) && (outflux < small_vel)) {
293  // do nothing
294  } else {
295  const Real alpha_fcf = influx/outflux; // flux correction factor
296 
297  correct_bdy_outflow_on_face(bdy_data_xlo, 0, true , domain, alpha_fcf, n);
298  correct_bdy_outflow_on_face(bdy_data_xhi, 0, false, domain, alpha_fcf, n);
299  correct_bdy_outflow_on_face(bdy_data_ylo, 1, true , domain, alpha_fcf, n);
300  correct_bdy_outflow_on_face(bdy_data_yhi, 1, false, domain, alpha_fcf, n);
301 
302  // For diagnostic purposes.
303  // Real influx_dbg = zero, outflux_dbg = zero;
304  // compute_influx_outflux_bdy(bdy_data_xlo, bdy_data_xhi,
305  // bdy_data_ylo, bdy_data_yhi,
306  // area_vec, geom, influx_dbg, outflux_dbg, n);
307  // amrex::Print() << " TOTAL BDY INFLUX / OUTFLOW AT WIDTH " << n << " " << influx_dbg << " " << outflux_dbg << "\n";
308  }
309  }
310 
311  multiply_by_rho0 = false;
312  scale_bdy_normal_by_rho0(rho0, bdy_data_xlo, 0, domain, domain_bcs_type_h, multiply_by_rho0);
313  scale_bdy_normal_by_rho0(rho0, bdy_data_xhi, 0, domain, domain_bcs_type_h, multiply_by_rho0);
314  scale_bdy_normal_by_rho0(rho0, bdy_data_ylo, 1, domain, domain_bcs_type_h, multiply_by_rho0);
315  scale_bdy_normal_by_rho0(rho0, bdy_data_yhi, 1, domain, domain_bcs_type_h, multiply_by_rho0);
316 }
compute_influx_outflux_bdy
void compute_influx_outflux_bdy(FArrayBox &bdy_data_xlo, FArrayBox &bdy_data_xhi, FArrayBox &bdy_data_ylo, FArrayBox &bdy_data_yhi, Array< MultiFab *, AMREX_SPACEDIM > &area_vec, const Geometry &geom, Real &influx, Real &outflux, const int n)
Definition: ERF_EnforceConstraintOnBdy.cpp:152
Here is the call graph for this function: