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 
)
138 {
139  BL_PROFILE_VAR("compute_influx_outflux_bdy()", computeInfluxOutfluxBdy);
140 
141  influx = zero, outflux = zero;
142 
143  const Box domain = geom.Domain();
144  const auto& domlo = lbound(domain);
145  const auto& domhi = ubound(domain);
146 
147  // Normal face area (of undistorted mesh)
148  const Real* a_dx = geom.CellSize();
149  const Real ds_x = a_dx[1];
150  const Real ds_y = a_dx[0];
151 
152  IntVect ngrow = {0,0,0};
153 
154  // X-dir
155  const auto& bdatxlo = bdy_data_xlo.const_array();
156  const auto& bdatxhi = bdy_data_xhi.const_array();
157  auto const& area_x = area_vec[0]->const_arrays();
158  influx += ds_x *
159  ParReduce(TypeList<ReduceOpSum>{},
160  TypeList<Real>{},
161  *area_vec[0], ngrow,
162  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
163  noexcept -> GpuTuple<Real>
164  {
165  if ( (i == domlo.x + n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxlo(i,j,k) > zero) {
166  return { std::abs(bdatxlo(i,j,k)) * area_x[box_no](i,j,k) };
167  } else if ( (i == domhi.x+1 - n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxhi(i,j,k) < zero) {
168  return { std::abs(bdatxhi(i,j,k)) * area_x[box_no](i,j,k) };
169  } else {
170  return { zero };
171  }
172  });
173 
174  outflux += ds_x *
175  ParReduce(TypeList<ReduceOpSum>{},
176  TypeList<Real>{},
177  *area_vec[0], ngrow,
178  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
179  noexcept -> GpuTuple<Real>
180  {
181  if (i == (domlo.x + n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxlo(i,j,k) < zero) {
182  return { std::abs(bdatxlo(i,j,k)) * area_x[box_no](i,j,k) };
183  } else if ( (i == domhi.x+1 - n) && (j >= domlo.y+n && j<= domhi.y-n) && bdatxhi(i,j,k) > zero) {
184  return { std::abs(bdatxhi(i,j,k)) * area_x[box_no](i,j,k) };
185  } else {
186  return { zero };
187  }
188  });
189 
190  // Y-dir
191  const auto& bdatylo = bdy_data_ylo.const_array();
192  const auto& bdatyhi = bdy_data_yhi.const_array();
193  auto const& area_y = area_vec[1]->const_arrays();
194  influx += ds_y *
195  ParReduce(TypeList<ReduceOpSum>{},
196  TypeList<Real>{},
197  *area_vec[1], ngrow,
198  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
199  noexcept -> GpuTuple<Real>
200  {
201  if ( (j == domlo.y + n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatylo(i,j,k) > zero) {
202  return { std::abs(bdatylo(i,j,k)) * area_y[box_no](i,j,k) };
203  } else if ( (j == domhi.y+1 - n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatyhi(i,j,k) < zero) {
204  return { std::abs(bdatyhi(i,j,k)) * area_y[box_no](i,j,k) };
205  } else {
206  return { zero };
207  }
208  });
209 
210  outflux += ds_y *
211  ParReduce(TypeList<ReduceOpSum>{},
212  TypeList<Real>{},
213  *area_vec[1], ngrow,
214  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
215  noexcept -> GpuTuple<Real>
216  {
217  if ( (j == domlo.y + n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatylo(i,j,k) < zero) {
218  return { std::abs(bdatylo(i,j,k)) * area_y[box_no](i,j,k) };
219  } else if ( (j == domhi.y+1 - n) && (i >= domlo.x+n && i<= domhi.x-n) && bdatyhi(i,j,k) > zero) {
220  return { std::abs(bdatyhi(i,j,k)) * area_y[box_no](i,j,k) };
221  } else {
222  return { zero };
223  }
224  });
225 
226  ParallelDescriptor::ReduceRealSum(influx);
227  ParallelDescriptor::ReduceRealSum(outflux);
228 }
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
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 
)
238 {
239  BL_PROFILE_VAR("enforceInOutSolvability_bdy()", enforceInOutSolvabilityBdy);
240 
241 #ifdef AMREX_USE_FLOAT
242  Real small_vel = Real(1.e-6);
243 #else
244  Real small_vel = Real(1.e-8);
245 #endif
246 
247  Box domain(geom.Domain());
248 
249  bool multiply_by_rho0 = true;
250  scale_bdy_normal_by_rho0(rho0, bdy_data_xlo, 0, domain, domain_bcs_type_h, multiply_by_rho0);
251  scale_bdy_normal_by_rho0(rho0, bdy_data_xhi, 0, domain, domain_bcs_type_h, multiply_by_rho0);
252  scale_bdy_normal_by_rho0(rho0, bdy_data_ylo, 1, domain, domain_bcs_type_h, multiply_by_rho0);
253  scale_bdy_normal_by_rho0(rho0, bdy_data_yhi, 1, domain, domain_bcs_type_h, multiply_by_rho0);
254 
255  int width =bdy_data_xlo.box().length(0);
256  AMREX_ASSERT(width == bdy_data_xhi.box().length(0));
257  AMREX_ASSERT(width == bdy_data_ylo.box().length(1));
258  AMREX_ASSERT(width == bdy_data_yhi.box().length(1));
259  // amrex::Print() << "width in enforce solvability " << width << std::endl;
260 
261  for (int n = 0; n < width; n++)
262  {
263  Real influx = zero, outflux = zero;
264  compute_influx_outflux_bdy(bdy_data_xlo, bdy_data_xhi,
265  bdy_data_ylo, bdy_data_yhi,
266  area_vec, geom, influx, outflux, n);
267 
268  // amrex::Print() << " TOTAL BDY INFLUX / OUTFLOW AT WIDTH " << n << " " << influx << " " << outflux << "\n";
269 
270  if ((influx > small_vel) && (outflux < small_vel)) {
271  Abort("Cannot enforce solvability on bdy_data, no outflow from the direction dependent boundaries");
272  } else if ((influx < small_vel) && (outflux < small_vel)) {
273  // do nothing
274  } else {
275  const Real alpha_fcf = influx/outflux; // flux correction factor
276 
277  correct_bdy_outflow_on_face(bdy_data_xlo, 0, true , domain, alpha_fcf, n);
278  correct_bdy_outflow_on_face(bdy_data_xhi, 0, false, domain, alpha_fcf, n);
279  correct_bdy_outflow_on_face(bdy_data_ylo, 1, true , domain, alpha_fcf, n);
280  correct_bdy_outflow_on_face(bdy_data_yhi, 1, false, domain, alpha_fcf, n);
281 
282  // For diagnostic purposes.
283  // Real influx_dbg = zero, outflux_dbg = zero;
284  // compute_influx_outflux_bdy(bdy_data_xlo, bdy_data_xhi,
285  // bdy_data_ylo, bdy_data_yhi,
286  // area_vec, geom, influx_dbg, outflux_dbg, n);
287  // amrex::Print() << " TOTAL BDY INFLUX / OUTFLOW AT WIDTH " << n << " " << influx_dbg << " " << outflux_dbg << "\n";
288  }
289  }
290 
291  multiply_by_rho0 = false;
292  scale_bdy_normal_by_rho0(rho0, bdy_data_xlo, 0, domain, domain_bcs_type_h, multiply_by_rho0);
293  scale_bdy_normal_by_rho0(rho0, bdy_data_xhi, 0, domain, domain_bcs_type_h, multiply_by_rho0);
294  scale_bdy_normal_by_rho0(rho0, bdy_data_ylo, 1, domain, domain_bcs_type_h, multiply_by_rho0);
295  scale_bdy_normal_by_rho0(rho0, bdy_data_yhi, 1, domain, domain_bcs_type_h, multiply_by_rho0);
296 }
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:132
Here is the call graph for this function: