From d3e9ac6a1f1477d3b27160d2cd2b1970eb547a9a Mon Sep 17 00:00:00 2001 From: asalmgren Date: Fri, 10 Jan 2025 02:00:13 +0000 Subject: [PATCH] Deployed from erf-model/ERF --- ERF_8H_source.html | 4 +- classERF.html | 1492 ++++++++++++++++++++++---------------------- 2 files changed, 750 insertions(+), 746 deletions(-) diff --git a/ERF_8H_source.html b/ERF_8H_source.html index f546f8eaa..151ef7512 100644 --- a/ERF_8H_source.html +++ b/ERF_8H_source.html @@ -1614,7 +1614,7 @@
ERF::init_sounding_ideal
static bool init_sounding_ideal
Definition: ERF.H:1073
ERF::ay_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > > ay_new
Definition: ERF.H:849
ERF::ERF
ERF()
Definition: ERF.cpp:87
-
ERF::init_Dirichlet_bc_data
void init_Dirichlet_bc_data(const std::string input_file)
Definition: ERF_InitBCs.cpp:599
+
ERF::init_Dirichlet_bc_data
void init_Dirichlet_bc_data(const std::string input_file)
Definition: ERF_InitBCs.cpp:600
ERF::~ERF
~ERF() override
ERF::restart_type
std::string restart_type
Definition: ERF.H:958
ERF::z_phys_nd_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > > z_phys_nd_src
Definition: ERF.H:840
@@ -1802,7 +1802,7 @@
ERF::SamplePointLogName
std::string SamplePointLogName(int i) const noexcept
The filename of the ith sampleptlog file.
Definition: ERF.H:1388
ERF::turbPert_update
void turbPert_update(const int lev, const amrex::Real dt)
Definition: ERF_InitTurbPert.cpp:12
ERF::InitData_post
void InitData_post()
Definition: ERF.cpp:654
-
ERF::refinement_criteria_setup
void refinement_criteria_setup()
Definition: ERF_Tagging.cpp:105
+
ERF::refinement_criteria_setup
void refinement_criteria_setup()
Definition: ERF_Tagging.cpp:108
ERF::nghost_eb_volume
static int nghost_eb_volume()
Definition: ERF.H:1407
ERF::metgrid_debug_dry
bool metgrid_debug_dry
Definition: ERF.H:1078
ERF::bndry_output_planes_interval
static int bndry_output_planes_interval
Definition: ERF.H:1100
diff --git a/classERF.html b/classERF.html index 9217ced89..5b40e03bd 100644 --- a/classERF.html +++ b/classERF.html @@ -5056,7 +5056,7 @@

ERF::yflux_imask
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > yflux_imask
Definition: ERF.H:898
ERF::az_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > > az_new
Definition: ERF.H:850
ERF::lsm_flux
amrex::Vector< amrex::Vector< amrex::MultiFab * > > lsm_flux
Definition: ERF.H:789
-
ERF::refinement_criteria_setup
void refinement_criteria_setup()
Definition: ERF_Tagging.cpp:105
+
ERF::refinement_criteria_setup
void refinement_criteria_setup()
Definition: ERF_Tagging.cpp:108
ERF::ax_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > > ax_src
Definition: ERF.H:842
ERF::zlevels_stag
amrex::Vector< amrex::Vector< amrex::Real > > zlevels_stag
Definition: ERF.H:829
ERF::lsm_data
amrex::Vector< amrex::Vector< amrex::MultiFab * > > lsm_data
Definition: ERF.H:788
@@ -5153,83 +5153,86 @@

19 
20  for (int j=0; j < ref_tags.size(); ++j)
21  {
-
22  std::unique_ptr<MultiFab> mf = std::make_unique<MultiFab>(grids[levc], dmap[levc], 1, 0);
-
23 
-
24  // This allows dynamic refinement based on the value of the density
-
25  if (ref_tags[j].Field() == "density")
-
26  {
-
27  MultiFab::Copy(*mf,vars_new[levc][Vars::cons],Rho_comp,0,1,0);
-
28 
-
29  // This allows dynamic refinement based on the value of qv
-
30  } else if ( ref_tags[j].Field() == "qv" ) {
-
31  MultiFab::Copy( *mf, vars_new[levc][Vars::cons], RhoQ1_comp, 0, 1, 0);
-
32  MultiFab::Divide(*mf, vars_new[levc][Vars::cons], Rho_comp, 0, 1, 0);
-
33 
-
34 
-
35  // This allows dynamic refinement based on the value of qc
-
36  } else if (ref_tags[j].Field() == "qc" ) {
-
37  MultiFab::Copy( *mf, vars_new[levc][Vars::cons], RhoQ2_comp, 0, 1, 0);
-
38  MultiFab::Divide(*mf, vars_new[levc][Vars::cons], Rho_comp, 0, 1, 0);
-
39 
-
40 
-
41  // This allows dynamic refinement based on the value of the scalar/pressure/theta
-
42  } else if ( (ref_tags[j].Field() == "scalar" ) ||
-
43  (ref_tags[j].Field() == "pressure") ||
-
44  (ref_tags[j].Field() == "theta" ) )
-
45  {
-
46  for (MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
-
47  {
-
48  const Box& bx = mfi.tilebox();
-
49  auto& dfab = (*mf)[mfi];
-
50  auto& sfab = vars_new[levc][Vars::cons][mfi];
-
51  if (ref_tags[j].Field() == "scalar") {
-
52  derived::erf_derscalar(bx, dfab, 0, 1, sfab, Geom(levc), time, nullptr, levc);
-
53  } else if (ref_tags[j].Field() == "theta") {
-
54  derived::erf_dertheta(bx, dfab, 0, 1, sfab, Geom(levc), time, nullptr, levc);
-
55  }
-
56  } // mfi
-
57 #ifdef ERF_USE_PARTICLES
-
58  } else {
-
59  //
-
60  // This allows dynamic refinement based on the number of particles per cell
-
61  //
-
62  // Note that we must count all the particles in levels both at and above the current,
-
63  // since otherwise, e.g., if the particles are all at level 1, counting particles at
-
64  // level 0 will not trigger refinement when regridding so level 1 will disappear,
-
65  // then come back at the next regridding
-
66  //
-
67  const auto& particles_namelist( particleData.getNames() );
-
68  mf->setVal(0.0);
-
69  for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++)
-
70  {
-
71  std::string tmp_string(particles_namelist[i]+"_count");
-
72  IntVect rr = IntVect::TheUnitVector();
-
73  if (ref_tags[j].Field() == tmp_string) {
-
74  for (int lev = levc; lev <= finest_level; lev++)
-
75  {
-
76  MultiFab temp_dat(grids[lev], dmap[lev], 1, 0); temp_dat.setVal(0);
-
77  particleData[particles_namelist[i]]->IncrementWithTotal(temp_dat, lev);
-
78 
-
79  MultiFab temp_dat_crse(grids[levc], dmap[levc], 1, 0); temp_dat_crse.setVal(0);
-
80 
-
81  if (lev == levc) {
-
82  MultiFab::Copy(*mf, temp_dat, 0, 0, 1, 0);
-
83  } else {
-
84  for (int d = 0; d < AMREX_SPACEDIM; d++) {
-
85  rr[d] *= ref_ratio[levc][d];
-
86  }
-
87  average_down(temp_dat, temp_dat_crse, 0, 1, rr);
-
88  MultiFab::Add(*mf, temp_dat_crse, 0, 0, 1, 0);
-
89  }
-
90  }
-
91  }
-
92  }
-
93 #endif
-
94  }
-
95 
-
96  ref_tags[j](tags,mf.get(),clearval,tagval,time,levc,geom[levc]);
-
97  } // loop over j
-
98 }
+
22  //
+
23  // This mf must have ghost cells because we may take differences between adjacent values
+
24  //
+
25  std::unique_ptr<MultiFab> mf = std::make_unique<MultiFab>(grids[levc], dmap[levc], 1, 1);
+
26 
+
27  // This allows dynamic refinement based on the value of the density
+
28  if (ref_tags[j].Field() == "density")
+
29  {
+
30  MultiFab::Copy(*mf,vars_new[levc][Vars::cons],Rho_comp,0,1,0);
+
31 
+
32  // This allows dynamic refinement based on the value of qv
+
33  } else if ( ref_tags[j].Field() == "qv" ) {
+
34  MultiFab::Copy( *mf, vars_new[levc][Vars::cons], RhoQ1_comp, 0, 1, 0);
+
35  MultiFab::Divide(*mf, vars_new[levc][Vars::cons], Rho_comp, 0, 1, 0);
+
36 
+
37 
+
38  // This allows dynamic refinement based on the value of qc
+
39  } else if (ref_tags[j].Field() == "qc" ) {
+
40  MultiFab::Copy( *mf, vars_new[levc][Vars::cons], RhoQ2_comp, 0, 1, 0);
+
41  MultiFab::Divide(*mf, vars_new[levc][Vars::cons], Rho_comp, 0, 1, 0);
+
42 
+
43 
+
44  // This allows dynamic refinement based on the value of the scalar/pressure/theta
+
45  } else if ( (ref_tags[j].Field() == "scalar" ) ||
+
46  (ref_tags[j].Field() == "pressure") ||
+
47  (ref_tags[j].Field() == "theta" ) )
+
48  {
+
49  for (MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+
50  {
+
51  const Box& bx = mfi.growntilebox();
+
52  auto& dfab = (*mf)[mfi];
+
53  auto& sfab = vars_new[levc][Vars::cons][mfi];
+
54  if (ref_tags[j].Field() == "scalar") {
+
55  derived::erf_derscalar(bx, dfab, 0, 1, sfab, Geom(levc), time, nullptr, levc);
+
56  } else if (ref_tags[j].Field() == "theta") {
+
57  derived::erf_dertheta(bx, dfab, 0, 1, sfab, Geom(levc), time, nullptr, levc);
+
58  }
+
59  } // mfi
+
60 #ifdef ERF_USE_PARTICLES
+
61  } else {
+
62  //
+
63  // This allows dynamic refinement based on the number of particles per cell
+
64  //
+
65  // Note that we must count all the particles in levels both at and above the current,
+
66  // since otherwise, e.g., if the particles are all at level 1, counting particles at
+
67  // level 0 will not trigger refinement when regridding so level 1 will disappear,
+
68  // then come back at the next regridding
+
69  //
+
70  const auto& particles_namelist( particleData.getNames() );
+
71  mf->setVal(0.0);
+
72  for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++)
+
73  {
+
74  std::string tmp_string(particles_namelist[i]+"_count");
+
75  IntVect rr = IntVect::TheUnitVector();
+
76  if (ref_tags[j].Field() == tmp_string) {
+
77  for (int lev = levc; lev <= finest_level; lev++)
+
78  {
+
79  MultiFab temp_dat(grids[lev], dmap[lev], 1, 0); temp_dat.setVal(0);
+
80  particleData[particles_namelist[i]]->IncrementWithTotal(temp_dat, lev);
+
81 
+
82  MultiFab temp_dat_crse(grids[levc], dmap[levc], 1, 0); temp_dat_crse.setVal(0);
+
83 
+
84  if (lev == levc) {
+
85  MultiFab::Copy(*mf, temp_dat, 0, 0, 1, 0);
+
86  } else {
+
87  for (int d = 0; d < AMREX_SPACEDIM; d++) {
+
88  rr[d] *= ref_ratio[levc][d];
+
89  }
+
90  average_down(temp_dat, temp_dat_crse, 0, 1, rr);
+
91  MultiFab::Add(*mf, temp_dat_crse, 0, 0, 1, 0);
+
92  }
+
93  }
+
94  }
+
95  }
+
96 #endif
+
97  }
+
98 
+
99  ref_tags[j](tags,mf.get(),clearval,tagval,time,levc,geom[levc]);
+
100  } // loop over j
+
101 }
ERF::ref_tags
static amrex::Vector< amrex::AMRErrorTag > ref_tags
Definition: ERF.H:1170
derived::erf_derscalar
void erf_derscalar(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: ERF_Derive.cpp:165
derived::erf_dertheta
void erf_dertheta(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: ERF_Derive.cpp:144
@@ -7114,379 +7117,380 @@

222  domain_bc_type[ori] = "SlipWall";
223 
224  Real rho_in;
-
225  if (pp.query("density", rho_in))
-
226  {
-
227  m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] = rho_in;
-
228  }
-
229 
-
230  Real theta_in;
-
231  if (pp.query("theta", theta_in))
-
232  {
-
233  m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] = theta_in*m_bc_extdir_vals[BCVars::Rho_bc_comp][ori];
-
234  }
-
235 
-
236  Real rho_grad_in;
-
237  if (pp.query("density_grad", rho_grad_in))
-
238  {
-
239  m_bc_neumann_vals[BCVars::Rho_bc_comp][ori] = rho_grad_in;
-
240  }
-
241 
-
242  Real theta_grad_in;
-
243  if (pp.query("theta_grad", theta_grad_in))
-
244  {
-
245  m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori] = theta_grad_in;
-
246  }
-
247  }
-
248  else if (bc_type == "most")
-
249  {
-
250  phys_bc_type[ori] = ERF_BC::MOST;
-
251  domain_bc_type[ori] = "MOST";
-
252  }
-
253  else
-
254  {
-
255  phys_bc_type[ori] = ERF_BC::undefined;
-
256  }
-
257 
-
258  if (geom[0].isPeriodic(ori.coordDir())) {
-
259  domain_bc_type[ori] = "Periodic";
-
260  if (phys_bc_type[ori] == ERF_BC::undefined)
-
261  {
-
262  phys_bc_type[ori] = ERF_BC::periodic;
-
263  } else {
-
264  Abort("Wrong BC type for periodic boundary");
-
265  }
-
266  }
-
267 
-
268  if (phys_bc_type[ori] == ERF_BC::undefined)
-
269  {
-
270  Print() << "BC Type specified for face " << bcid << " is " << bc_type_in << std::endl;
-
271  Abort("This BC type is unknown");
-
272  }
-
273  };
-
274 
-
275  f("xlo", Orientation(Direction::x,Orientation::low));
-
276  f("xhi", Orientation(Direction::x,Orientation::high));
-
277  f("ylo", Orientation(Direction::y,Orientation::low));
-
278  f("yhi", Orientation(Direction::y,Orientation::high));
-
279  f("zlo", Orientation(Direction::z,Orientation::low));
-
280  f("zhi", Orientation(Direction::z,Orientation::high));
-
281 
-
282  // *****************************************************************************
-
283  //
-
284  // Here we translate the physical boundary conditions -- one type per face --
-
285  // into logical boundary conditions for each velocity component
-
286  //
-
287  // *****************************************************************************
-
288  {
-
289  domain_bcs_type.resize(AMREX_SPACEDIM+NBCVAR_max);
-
290  domain_bcs_type_d.resize(AMREX_SPACEDIM+NBCVAR_max);
-
291 
-
292  for (OrientationIter oit; oit; ++oit) {
-
293  Orientation ori = oit();
-
294  int dir = ori.coordDir();
-
295  Orientation::Side side = ori.faceDir();
-
296  auto const bct = phys_bc_type[ori];
-
297  if ( bct == ERF_BC::symmetry )
-
298  {
-
299  if (side == Orientation::low) {
-
300  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
301  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::reflect_even);
-
302  }
-
303  domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::reflect_odd);
-
304  } else {
-
305  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
306  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::reflect_even);
-
307  }
-
308  domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::reflect_odd);
-
309  }
-
310  }
-
311  else if (bct == ERF_BC::outflow or bct == ERF_BC::ho_outflow )
-
312  {
-
313  if (side == Orientation::low) {
-
314  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
315  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::foextrap);
-
316  }
-
317  domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::neumann_int);
-
318  } else {
-
319  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
320  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::foextrap);
-
321  }
-
322  domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::neumann_int);
-
323  }
-
324  }
-
325  else if (bct == ERF_BC::open)
-
326  {
-
327  if (side == Orientation::low) {
-
328  for (int i = 0; i < AMREX_SPACEDIM; i++)
-
329  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::open);
-
330  } else {
-
331  for (int i = 0; i < AMREX_SPACEDIM; i++)
-
332  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::open);
-
333  }
-
334  }
-
335  else if (bct == ERF_BC::inflow)
-
336  {
-
337  if (side == Orientation::low) {
-
338  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
339  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir);
-
340  if (input_bndry_planes && dir < 2 && m_r2d->ingested_velocity()) {
-
341  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir_ingested);
-
342  }
-
343  }
-
344  } else {
-
345  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
346  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir);
-
347  if (input_bndry_planes && dir < 2 && m_r2d->ingested_velocity()) {
-
348  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir_ingested);
-
349  }
-
350  }
-
351  }
-
352  }
-
353  else if (bct == ERF_BC::no_slip_wall)
-
354  {
-
355  if (side == Orientation::low) {
-
356  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
357  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir);
-
358  }
-
359  } else {
-
360  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
361  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir);
-
362  }
-
363  }
-
364  }
-
365  else if (bct == ERF_BC::slip_wall)
-
366  {
-
367  if (side == Orientation::low) {
-
368  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
369  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::foextrap);
-
370  }
-
371  // Only normal direction has ext_dir
-
372  domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::ext_dir);
-
373 
-
374  } else {
-
375  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
376  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::foextrap);
-
377  }
-
378  // Only normal direction has ext_dir
-
379  domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::ext_dir);
-
380  }
-
381  }
-
382  else if (bct == ERF_BC::periodic)
-
383  {
-
384  if (side == Orientation::low) {
-
385  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
386  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::int_dir);
-
387  }
-
388  } else {
-
389  for (int i = 0; i < AMREX_SPACEDIM; i++) {
-
390  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::int_dir);
-
391  }
-
392  }
-
393  }
-
394  else if ( bct == ERF_BC::MOST )
-
395  {
-
396  AMREX_ALWAYS_ASSERT(dir == 2 && side == Orientation::low);
-
397  domain_bcs_type[BCVars::xvel_bc+0].setLo(dir, ERFBCType::reflect_odd);
-
398  domain_bcs_type[BCVars::xvel_bc+1].setLo(dir, ERFBCType::reflect_odd);
-
399  domain_bcs_type[BCVars::xvel_bc+2].setLo(dir, ERFBCType::ext_dir);
-
400  }
-
401  }
-
402  }
-
403 
-
404  // *****************************************************************************
-
405  //
-
406  // Here we translate the physical boundary conditions -- one type per face --
-
407  // into logical boundary conditions for each cell-centered variable
-
408  // (including the base state variables)
-
409  // NOTE: all "scalars" share the same type of boundary condition
-
410  //
-
411  // *****************************************************************************
-
412  {
-
413  for (OrientationIter oit; oit; ++oit) {
-
414  Orientation ori = oit();
-
415  int dir = ori.coordDir();
-
416  Orientation::Side side = ori.faceDir();
-
417  auto const bct = phys_bc_type[ori];
-
418  if ( bct == ERF_BC::symmetry )
-
419  {
-
420  if (side == Orientation::low) {
-
421  for (int i = 0; i < NBCVAR_max; i++) {
-
422  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::reflect_even);
-
423  }
-
424  } else {
-
425  for (int i = 0; i < NBCVAR_max; i++) {
-
426  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::reflect_even);
-
427  }
-
428  }
-
429  }
-
430  else if ( bct == ERF_BC::outflow )
-
431  {
-
432  if (side == Orientation::low) {
-
433  for (int i = 0; i < NBCVAR_max; i++) {
-
434  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
-
435  }
-
436  } else {
-
437  for (int i = 0; i < NBCVAR_max; i++) {
-
438  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
-
439  }
-
440  }
-
441  }
-
442  else if ( bct == ERF_BC::ho_outflow )
-
443  {
-
444  if (side == Orientation::low) {
-
445  for (int i = 0; i < NBCVAR_max; i++) {
-
446  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::hoextrapcc);
-
447  }
-
448  } else {
-
449  for (int i = 0; i < NBCVAR_max; i++) {
-
450  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::hoextrapcc);
-
451  }
-
452  }
-
453  }
-
454  else if ( bct == ERF_BC::open )
-
455  {
-
456  if (side == Orientation::low) {
-
457  for (int i = 0; i < NBCVAR_max; i++)
-
458  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::open);
-
459  } else {
-
460  for (int i = 0; i < NBCVAR_max; i++)
-
461  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::open);
-
462  }
-
463  }
-
464  else if ( bct == ERF_BC::no_slip_wall )
-
465  {
-
466  if (side == Orientation::low) {
-
467  for (int i = 0; i < NBCVAR_max; i++) {
-
468  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
-
469  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
-
470  if (rho_read) {
-
471  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir);
-
472  } else {
-
473  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir_prim);
-
474  }
-
475  }
-
476  }
-
477  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
-
478  domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::neumann);
-
479  }
-
480  } else {
-
481  for (int i = 0; i < NBCVAR_max; i++) {
-
482  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
-
483  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
-
484  if (rho_read) {
-
485  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir);
-
486  } else {
-
487  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir_prim);
-
488  }
-
489  }
-
490  }
-
491  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
-
492  domain_bcs_type[BCVars::RhoTheta_bc_comp].setHi(dir, ERFBCType::neumann);
-
493  }
-
494  }
-
495  }
-
496  else if (bct == ERF_BC::slip_wall)
-
497  {
-
498  if (side == Orientation::low) {
-
499  for (int i = 0; i < NBCVAR_max; i++) {
-
500  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
-
501  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
-
502  if (rho_read) {
-
503  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir);
-
504  } else {
-
505  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir_prim);
-
506  }
-
507  }
-
508  }
-
509  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
-
510  domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::neumann);
-
511  }
-
512  if (std::abs(m_bc_neumann_vals[BCVars::Rho_bc_comp][ori]) > 0.) {
-
513  domain_bcs_type[BCVars::Rho_bc_comp].setLo(dir, ERFBCType::neumann);
-
514  }
-
515  } else {
-
516  for (int i = 0; i < NBCVAR_max; i++) {
-
517  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
-
518  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
-
519  if (rho_read) {
-
520  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir);
-
521  } else {
-
522  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir_prim);
-
523  }
-
524  }
-
525  }
-
526  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
-
527  domain_bcs_type[BCVars::RhoTheta_bc_comp].setHi(dir, ERFBCType::neumann);
-
528  }
-
529  if (std::abs(m_bc_neumann_vals[BCVars::Rho_bc_comp][ori]) > 0.) {
-
530  domain_bcs_type[BCVars::Rho_bc_comp].setHi(dir, ERFBCType::neumann);
-
531  }
-
532  }
-
533  }
-
534  else if (bct == ERF_BC::inflow)
-
535  {
-
536  if (side == Orientation::low) {
-
537  for (int i = 0; i < NBCVAR_max; i++) {
-
538  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir);
-
539  if (input_bndry_planes && dir < 2 && (
-
540  ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) ||
-
541  ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) ||
-
542  ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) ||
-
543  ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) ||
-
544  ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) ||
-
545  ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() )) )
-
546  {
-
547  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir_ingested);
-
548  }
-
549  else if (m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] == 0) {
-
550  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
-
551  }
-
552  }
-
553  } else {
-
554  for (int i = 0; i < NBCVAR_max; i++) {
-
555  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir);
-
556  if (input_bndry_planes && dir < 2 && (
-
557  ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) ||
-
558  ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) ||
-
559  ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) ||
-
560  ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) ||
-
561  ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) ||
-
562  ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() )
-
563  ) )
-
564  {
-
565  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir_ingested);
-
566  }
-
567  else if (m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] == 0) {
-
568  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
-
569  }
-
570  }
-
571  }
-
572  }
-
573  else if (bct == ERF_BC::periodic)
-
574  {
-
575  if (side == Orientation::low) {
-
576  for (int i = 0; i < NBCVAR_max; i++) {
-
577  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::int_dir);
-
578  }
-
579  } else {
-
580  for (int i = 0; i < NBCVAR_max; i++) {
-
581  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::int_dir);
-
582  }
-
583  }
-
584  }
-
585  else if ( bct == ERF_BC::MOST )
-
586  {
-
587  AMREX_ALWAYS_ASSERT(dir == 2 && side == Orientation::low);
-
588  for (int i = 0; i < NBCVAR_max; i++) {
-
589  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
-
590  }
-
591  }
-
592  }
-
593  }
-
594 
-
595  // NOTE: Gpu:copy is a wrapper to htod_memcpy (GPU) or memcpy (CPU) and is a blocking comm
-
596  Gpu::copy(Gpu::hostToDevice, domain_bcs_type.begin(), domain_bcs_type.end(), domain_bcs_type_d.begin());
-
597 }
+
225  rho_read = pp.query("density", rho_in);
+
226  if (rho_read)
+
227  {
+
228  m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] = rho_in;
+
229  }
+
230 
+
231  Real theta_in;
+
232  if (pp.query("theta", theta_in))
+
233  {
+
234  m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] = theta_in*m_bc_extdir_vals[BCVars::Rho_bc_comp][ori];
+
235  }
+
236 
+
237  Real rho_grad_in;
+
238  if (pp.query("density_grad", rho_grad_in))
+
239  {
+
240  m_bc_neumann_vals[BCVars::Rho_bc_comp][ori] = rho_grad_in;
+
241  }
+
242 
+
243  Real theta_grad_in;
+
244  if (pp.query("theta_grad", theta_grad_in))
+
245  {
+
246  m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori] = theta_grad_in;
+
247  }
+
248  }
+
249  else if (bc_type == "most")
+
250  {
+
251  phys_bc_type[ori] = ERF_BC::MOST;
+
252  domain_bc_type[ori] = "MOST";
+
253  }
+
254  else
+
255  {
+
256  phys_bc_type[ori] = ERF_BC::undefined;
+
257  }
+
258 
+
259  if (geom[0].isPeriodic(ori.coordDir())) {
+
260  domain_bc_type[ori] = "Periodic";
+
261  if (phys_bc_type[ori] == ERF_BC::undefined)
+
262  {
+
263  phys_bc_type[ori] = ERF_BC::periodic;
+
264  } else {
+
265  Abort("Wrong BC type for periodic boundary");
+
266  }
+
267  }
+
268 
+
269  if (phys_bc_type[ori] == ERF_BC::undefined)
+
270  {
+
271  Print() << "BC Type specified for face " << bcid << " is " << bc_type_in << std::endl;
+
272  Abort("This BC type is unknown");
+
273  }
+
274  };
+
275 
+
276  f("xlo", Orientation(Direction::x,Orientation::low));
+
277  f("xhi", Orientation(Direction::x,Orientation::high));
+
278  f("ylo", Orientation(Direction::y,Orientation::low));
+
279  f("yhi", Orientation(Direction::y,Orientation::high));
+
280  f("zlo", Orientation(Direction::z,Orientation::low));
+
281  f("zhi", Orientation(Direction::z,Orientation::high));
+
282 
+
283  // *****************************************************************************
+
284  //
+
285  // Here we translate the physical boundary conditions -- one type per face --
+
286  // into logical boundary conditions for each velocity component
+
287  //
+
288  // *****************************************************************************
+
289  {
+
290  domain_bcs_type.resize(AMREX_SPACEDIM+NBCVAR_max);
+
291  domain_bcs_type_d.resize(AMREX_SPACEDIM+NBCVAR_max);
+
292 
+
293  for (OrientationIter oit; oit; ++oit) {
+
294  Orientation ori = oit();
+
295  int dir = ori.coordDir();
+
296  Orientation::Side side = ori.faceDir();
+
297  auto const bct = phys_bc_type[ori];
+
298  if ( bct == ERF_BC::symmetry )
+
299  {
+
300  if (side == Orientation::low) {
+
301  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
302  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::reflect_even);
+
303  }
+
304  domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::reflect_odd);
+
305  } else {
+
306  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
307  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::reflect_even);
+
308  }
+
309  domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::reflect_odd);
+
310  }
+
311  }
+
312  else if (bct == ERF_BC::outflow or bct == ERF_BC::ho_outflow )
+
313  {
+
314  if (side == Orientation::low) {
+
315  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
316  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::foextrap);
+
317  }
+
318  domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::neumann_int);
+
319  } else {
+
320  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
321  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::foextrap);
+
322  }
+
323  domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::neumann_int);
+
324  }
+
325  }
+
326  else if (bct == ERF_BC::open)
+
327  {
+
328  if (side == Orientation::low) {
+
329  for (int i = 0; i < AMREX_SPACEDIM; i++)
+
330  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::open);
+
331  } else {
+
332  for (int i = 0; i < AMREX_SPACEDIM; i++)
+
333  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::open);
+
334  }
+
335  }
+
336  else if (bct == ERF_BC::inflow)
+
337  {
+
338  if (side == Orientation::low) {
+
339  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
340  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir);
+
341  if (input_bndry_planes && dir < 2 && m_r2d->ingested_velocity()) {
+
342  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir_ingested);
+
343  }
+
344  }
+
345  } else {
+
346  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
347  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir);
+
348  if (input_bndry_planes && dir < 2 && m_r2d->ingested_velocity()) {
+
349  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir_ingested);
+
350  }
+
351  }
+
352  }
+
353  }
+
354  else if (bct == ERF_BC::no_slip_wall)
+
355  {
+
356  if (side == Orientation::low) {
+
357  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
358  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir);
+
359  }
+
360  } else {
+
361  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
362  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir);
+
363  }
+
364  }
+
365  }
+
366  else if (bct == ERF_BC::slip_wall)
+
367  {
+
368  if (side == Orientation::low) {
+
369  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
370  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::foextrap);
+
371  }
+
372  // Only normal direction has ext_dir
+
373  domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::ext_dir);
+
374 
+
375  } else {
+
376  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
377  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::foextrap);
+
378  }
+
379  // Only normal direction has ext_dir
+
380  domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::ext_dir);
+
381  }
+
382  }
+
383  else if (bct == ERF_BC::periodic)
+
384  {
+
385  if (side == Orientation::low) {
+
386  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
387  domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::int_dir);
+
388  }
+
389  } else {
+
390  for (int i = 0; i < AMREX_SPACEDIM; i++) {
+
391  domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::int_dir);
+
392  }
+
393  }
+
394  }
+
395  else if ( bct == ERF_BC::MOST )
+
396  {
+
397  AMREX_ALWAYS_ASSERT(dir == 2 && side == Orientation::low);
+
398  domain_bcs_type[BCVars::xvel_bc+0].setLo(dir, ERFBCType::reflect_odd);
+
399  domain_bcs_type[BCVars::xvel_bc+1].setLo(dir, ERFBCType::reflect_odd);
+
400  domain_bcs_type[BCVars::xvel_bc+2].setLo(dir, ERFBCType::ext_dir);
+
401  }
+
402  }
+
403  }
+
404 
+
405  // *****************************************************************************
+
406  //
+
407  // Here we translate the physical boundary conditions -- one type per face --
+
408  // into logical boundary conditions for each cell-centered variable
+
409  // (including the base state variables)
+
410  // NOTE: all "scalars" share the same type of boundary condition
+
411  //
+
412  // *****************************************************************************
+
413  {
+
414  for (OrientationIter oit; oit; ++oit) {
+
415  Orientation ori = oit();
+
416  int dir = ori.coordDir();
+
417  Orientation::Side side = ori.faceDir();
+
418  auto const bct = phys_bc_type[ori];
+
419  if ( bct == ERF_BC::symmetry )
+
420  {
+
421  if (side == Orientation::low) {
+
422  for (int i = 0; i < NBCVAR_max; i++) {
+
423  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::reflect_even);
+
424  }
+
425  } else {
+
426  for (int i = 0; i < NBCVAR_max; i++) {
+
427  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::reflect_even);
+
428  }
+
429  }
+
430  }
+
431  else if ( bct == ERF_BC::outflow )
+
432  {
+
433  if (side == Orientation::low) {
+
434  for (int i = 0; i < NBCVAR_max; i++) {
+
435  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
+
436  }
+
437  } else {
+
438  for (int i = 0; i < NBCVAR_max; i++) {
+
439  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
+
440  }
+
441  }
+
442  }
+
443  else if ( bct == ERF_BC::ho_outflow )
+
444  {
+
445  if (side == Orientation::low) {
+
446  for (int i = 0; i < NBCVAR_max; i++) {
+
447  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::hoextrapcc);
+
448  }
+
449  } else {
+
450  for (int i = 0; i < NBCVAR_max; i++) {
+
451  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::hoextrapcc);
+
452  }
+
453  }
+
454  }
+
455  else if ( bct == ERF_BC::open )
+
456  {
+
457  if (side == Orientation::low) {
+
458  for (int i = 0; i < NBCVAR_max; i++)
+
459  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::open);
+
460  } else {
+
461  for (int i = 0; i < NBCVAR_max; i++)
+
462  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::open);
+
463  }
+
464  }
+
465  else if ( bct == ERF_BC::no_slip_wall )
+
466  {
+
467  if (side == Orientation::low) {
+
468  for (int i = 0; i < NBCVAR_max; i++) {
+
469  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
+
470  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
+
471  if (rho_read) {
+
472  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir);
+
473  } else {
+
474  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir_prim);
+
475  }
+
476  }
+
477  }
+
478  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
+
479  domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::neumann);
+
480  }
+
481  } else {
+
482  for (int i = 0; i < NBCVAR_max; i++) {
+
483  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
+
484  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
+
485  if (rho_read) {
+
486  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir);
+
487  } else {
+
488  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir_prim);
+
489  }
+
490  }
+
491  }
+
492  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
+
493  domain_bcs_type[BCVars::RhoTheta_bc_comp].setHi(dir, ERFBCType::neumann);
+
494  }
+
495  }
+
496  }
+
497  else if (bct == ERF_BC::slip_wall)
+
498  {
+
499  if (side == Orientation::low) {
+
500  for (int i = 0; i < NBCVAR_max; i++) {
+
501  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
+
502  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
+
503  if (rho_read) {
+
504  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir);
+
505  } else {
+
506  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir_prim);
+
507  }
+
508  }
+
509  }
+
510  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
+
511  domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::neumann);
+
512  }
+
513  if (std::abs(m_bc_neumann_vals[BCVars::Rho_bc_comp][ori]) > 0.) {
+
514  domain_bcs_type[BCVars::Rho_bc_comp].setLo(dir, ERFBCType::neumann);
+
515  }
+
516  } else {
+
517  for (int i = 0; i < NBCVAR_max; i++) {
+
518  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
+
519  if (m_bc_extdir_vals[BCVars::cons_bc+i][ori] != cons_dir_init[BCVars::cons_bc+i]) {
+
520  if (rho_read) {
+
521  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir);
+
522  } else {
+
523  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir_prim);
+
524  }
+
525  }
+
526  }
+
527  if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) {
+
528  domain_bcs_type[BCVars::RhoTheta_bc_comp].setHi(dir, ERFBCType::neumann);
+
529  }
+
530  if (std::abs(m_bc_neumann_vals[BCVars::Rho_bc_comp][ori]) > 0.) {
+
531  domain_bcs_type[BCVars::Rho_bc_comp].setHi(dir, ERFBCType::neumann);
+
532  }
+
533  }
+
534  }
+
535  else if (bct == ERF_BC::inflow)
+
536  {
+
537  if (side == Orientation::low) {
+
538  for (int i = 0; i < NBCVAR_max; i++) {
+
539  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir);
+
540  if (input_bndry_planes && dir < 2 && (
+
541  ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) ||
+
542  ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) ||
+
543  ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) ||
+
544  ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) ||
+
545  ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) ||
+
546  ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() )) )
+
547  {
+
548  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir_ingested);
+
549  }
+
550  else if (m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] == 0) {
+
551  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
+
552  }
+
553  }
+
554  } else {
+
555  for (int i = 0; i < NBCVAR_max; i++) {
+
556  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir);
+
557  if (input_bndry_planes && dir < 2 && (
+
558  ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) ||
+
559  ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) ||
+
560  ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) ||
+
561  ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) ||
+
562  ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) ||
+
563  ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() )
+
564  ) )
+
565  {
+
566  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir_ingested);
+
567  }
+
568  else if (m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] == 0) {
+
569  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap);
+
570  }
+
571  }
+
572  }
+
573  }
+
574  else if (bct == ERF_BC::periodic)
+
575  {
+
576  if (side == Orientation::low) {
+
577  for (int i = 0; i < NBCVAR_max; i++) {
+
578  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::int_dir);
+
579  }
+
580  } else {
+
581  for (int i = 0; i < NBCVAR_max; i++) {
+
582  domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::int_dir);
+
583  }
+
584  }
+
585  }
+
586  else if ( bct == ERF_BC::MOST )
+
587  {
+
588  AMREX_ALWAYS_ASSERT(dir == 2 && side == Orientation::low);
+
589  for (int i = 0; i < NBCVAR_max; i++) {
+
590  domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap);
+
591  }
+
592  }
+
593  }
+
594  }
+
595 
+
596  // NOTE: Gpu:copy is a wrapper to htod_memcpy (GPU) or memcpy (CPU) and is a blocking comm
+
597  Gpu::copy(Gpu::hostToDevice, domain_bcs_type.begin(), domain_bcs_type.end(), domain_bcs_type_d.begin());
+
598 }
NBCVAR_max
#define NBCVAR_max
Definition: ERF_IndexDefines.H:29
ERF_BC::slip_wall
@ slip_wall
ERF_BC::symmetry
@ symmetry
@@ -7496,7 +7500,7 @@

ERF_BC::outflow
@ outflow
ERF_BC::inflow
@ inflow
ERF_BC::periodic
@ periodic
-
ERF::init_Dirichlet_bc_data
void init_Dirichlet_bc_data(const std::string input_file)
Definition: ERF_InitBCs.cpp:599
+
ERF::init_Dirichlet_bc_data
void init_Dirichlet_bc_data(const std::string input_file)
Definition: ERF_InitBCs.cpp:600
ERF::m_bc_neumann_vals
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_neumann_vals
Definition: ERF.H:891
BCVars::RhoQ6_bc_comp
@ RhoQ6_bc_comp
Definition: ERF_IndexDefines.H:86
BCVars::RhoQ1_bc_comp
@ RhoQ1_bc_comp
Definition: ERF_IndexDefines.H:81
@@ -7671,91 +7675,91 @@

-
600 {
-
601  const bool use_terrain = (solverChoice.terrain_type != TerrainType::None);
-
602 
-
603  // Read the dirichlet_input file
-
604  Print() << "dirichlet_input file location : " << input_file << std::endl;
-
605  std::ifstream input_reader(input_file);
-
606  if (!input_reader.is_open()) {
-
607  amrex::Abort("Error opening the dirichlet_input file.\n");
-
608  }
-
609 
-
610  Print() << "Successfully opened the dirichlet_input file. Now reading... " << std::endl;
-
611  std::string line;
-
612 
-
613  // Size of Ninp (number of z points in input file)
-
614  Vector<Real> z_inp_tmp, u_inp_tmp, v_inp_tmp, w_inp_tmp;
-
615 
-
616  const int klo = geom[0].Domain().smallEnd()[2];
-
617  const int khi = geom[0].Domain().bigEnd()[2];
-
618 
-
619  const Real zbot = (use_terrain) ? zlevels_stag[0][klo] : geom[0].ProbLo(2);
-
620  const Real ztop = (use_terrain) ? zlevels_stag[0][khi+1] : geom[0].ProbHi(2);
-
621 
-
622  // Add surface
-
623  z_inp_tmp.push_back(zbot); // height above sea level [m]
-
624  u_inp_tmp.push_back(0.);
-
625  v_inp_tmp.push_back(0.);
-
626  w_inp_tmp.push_back(0.);
-
627 
-
628  // Read the vertical profile at each given height
-
629  Real z, u, v, w;
-
630  while(std::getline(input_reader, line)) {
-
631  std::istringstream iss_z(line);
-
632  iss_z >> z >> u >> v >> w;
-
633  if (z == zbot) {
-
634  u_inp_tmp[0] = u;
-
635  v_inp_tmp[0] = v;
-
636  w_inp_tmp[0] = w;
-
637  } else {
-
638  AMREX_ALWAYS_ASSERT(z > z_inp_tmp[z_inp_tmp.size()-1]); // sounding is increasing in height
-
639  z_inp_tmp.push_back(z);
-
640  u_inp_tmp.push_back(u);
-
641  v_inp_tmp.push_back(v);
-
642  w_inp_tmp.push_back(w);
-
643  if (z >= ztop) break;
-
644  }
-
645  }
-
646 
-
647  amrex::Print() << "Successfully read and interpolated the dirichlet_input file..." << std::endl;
-
648  input_reader.close();
-
649 
-
650  for (int lev = 0; lev <= max_level; lev++) {
-
651 
-
652  const int Nz = geom[lev].Domain().size()[2];
-
653  const Real dz = geom[lev].CellSize()[2];
-
654 
-
655  // Size of Nz (domain grid)
-
656  Vector<Real> zcc_inp(Nz );
-
657  Vector<Real> znd_inp(Nz+1);
-
658  Vector<Real> u_inp(Nz ); xvel_bc_data[lev].resize(Nz ,0.0);
-
659  Vector<Real> v_inp(Nz ); yvel_bc_data[lev].resize(Nz ,0.0);
-
660  Vector<Real> w_inp(Nz+1); zvel_bc_data[lev].resize(Nz+1,0.0);
-
661 
-
662  // At this point, we have an input from zbot up to
-
663  // z_inp_tmp[N-1] >= ztop. Now, interpolate to grid level 0 heights
-
664  const int Ninp = z_inp_tmp.size();
-
665  for (int k(0); k<Nz; ++k) {
-
666  zcc_inp[k] = (use_terrain) ? 0.5 * (zlevels_stag[0][k] + zlevels_stag[0][k+1])
-
667  : zbot + (k + 0.5) * dz;
-
668  znd_inp[k] = (use_terrain) ? zlevels_stag[0][k+1] : zbot + (k) * dz;
-
669  u_inp[k] = interpolate_1d(z_inp_tmp.dataPtr(), u_inp_tmp.dataPtr(), zcc_inp[k], Ninp);
-
670  v_inp[k] = interpolate_1d(z_inp_tmp.dataPtr(), v_inp_tmp.dataPtr(), zcc_inp[k], Ninp);
-
671  w_inp[k] = interpolate_1d(z_inp_tmp.dataPtr(), w_inp_tmp.dataPtr(), znd_inp[k], Ninp);
-
672  }
-
673  znd_inp[Nz] = ztop;
-
674  w_inp[Nz] = interpolate_1d(z_inp_tmp.dataPtr(), w_inp_tmp.dataPtr(), ztop, Ninp);
-
675 
-
676  // Copy host data to the device
-
677  Gpu::copy(Gpu::hostToDevice, u_inp.begin(), u_inp.end(), xvel_bc_data[lev].begin());
-
678  Gpu::copy(Gpu::hostToDevice, v_inp.begin(), v_inp.end(), yvel_bc_data[lev].begin());
-
679  Gpu::copy(Gpu::hostToDevice, w_inp.begin(), w_inp.end(), zvel_bc_data[lev].begin());
-
680 
-
681  // NOTE: These device vectors are passed to the PhysBC constructors when that
-
682  // class is instantiated in ERF_MakeNewArrays.cpp.
-
683  } // lev
-
684 }
+
601 {
+
602  const bool use_terrain = (solverChoice.terrain_type != TerrainType::None);
+
603 
+
604  // Read the dirichlet_input file
+
605  Print() << "dirichlet_input file location : " << input_file << std::endl;
+
606  std::ifstream input_reader(input_file);
+
607  if (!input_reader.is_open()) {
+
608  amrex::Abort("Error opening the dirichlet_input file.\n");
+
609  }
+
610 
+
611  Print() << "Successfully opened the dirichlet_input file. Now reading... " << std::endl;
+
612  std::string line;
+
613 
+
614  // Size of Ninp (number of z points in input file)
+
615  Vector<Real> z_inp_tmp, u_inp_tmp, v_inp_tmp, w_inp_tmp;
+
616 
+
617  const int klo = geom[0].Domain().smallEnd()[2];
+
618  const int khi = geom[0].Domain().bigEnd()[2];
+
619 
+
620  const Real zbot = (use_terrain) ? zlevels_stag[0][klo] : geom[0].ProbLo(2);
+
621  const Real ztop = (use_terrain) ? zlevels_stag[0][khi+1] : geom[0].ProbHi(2);
+
622 
+
623  // Add surface
+
624  z_inp_tmp.push_back(zbot); // height above sea level [m]
+
625  u_inp_tmp.push_back(0.);
+
626  v_inp_tmp.push_back(0.);
+
627  w_inp_tmp.push_back(0.);
+
628 
+
629  // Read the vertical profile at each given height
+
630  Real z, u, v, w;
+
631  while(std::getline(input_reader, line)) {
+
632  std::istringstream iss_z(line);
+
633  iss_z >> z >> u >> v >> w;
+
634  if (z == zbot) {
+
635  u_inp_tmp[0] = u;
+
636  v_inp_tmp[0] = v;
+
637  w_inp_tmp[0] = w;
+
638  } else {
+
639  AMREX_ALWAYS_ASSERT(z > z_inp_tmp[z_inp_tmp.size()-1]); // sounding is increasing in height
+
640  z_inp_tmp.push_back(z);
+
641  u_inp_tmp.push_back(u);
+
642  v_inp_tmp.push_back(v);
+
643  w_inp_tmp.push_back(w);
+
644  if (z >= ztop) break;
+
645  }
+
646  }
+
647 
+
648  amrex::Print() << "Successfully read and interpolated the dirichlet_input file..." << std::endl;
+
649  input_reader.close();
+
650 
+
651  for (int lev = 0; lev <= max_level; lev++) {
+
652 
+
653  const int Nz = geom[lev].Domain().size()[2];
+
654  const Real dz = geom[lev].CellSize()[2];
+
655 
+
656  // Size of Nz (domain grid)
+
657  Vector<Real> zcc_inp(Nz );
+
658  Vector<Real> znd_inp(Nz+1);
+
659  Vector<Real> u_inp(Nz ); xvel_bc_data[lev].resize(Nz ,0.0);
+
660  Vector<Real> v_inp(Nz ); yvel_bc_data[lev].resize(Nz ,0.0);
+
661  Vector<Real> w_inp(Nz+1); zvel_bc_data[lev].resize(Nz+1,0.0);
+
662 
+
663  // At this point, we have an input from zbot up to
+
664  // z_inp_tmp[N-1] >= ztop. Now, interpolate to grid level 0 heights
+
665  const int Ninp = z_inp_tmp.size();
+
666  for (int k(0); k<Nz; ++k) {
+
667  zcc_inp[k] = (use_terrain) ? 0.5 * (zlevels_stag[0][k] + zlevels_stag[0][k+1])
+
668  : zbot + (k + 0.5) * dz;
+
669  znd_inp[k] = (use_terrain) ? zlevels_stag[0][k+1] : zbot + (k) * dz;
+
670  u_inp[k] = interpolate_1d(z_inp_tmp.dataPtr(), u_inp_tmp.dataPtr(), zcc_inp[k], Ninp);
+
671  v_inp[k] = interpolate_1d(z_inp_tmp.dataPtr(), v_inp_tmp.dataPtr(), zcc_inp[k], Ninp);
+
672  w_inp[k] = interpolate_1d(z_inp_tmp.dataPtr(), w_inp_tmp.dataPtr(), znd_inp[k], Ninp);
+
673  }
+
674  znd_inp[Nz] = ztop;
+
675  w_inp[Nz] = interpolate_1d(z_inp_tmp.dataPtr(), w_inp_tmp.dataPtr(), ztop, Ninp);
+
676 
+
677  // Copy host data to the device
+
678  Gpu::copy(Gpu::hostToDevice, u_inp.begin(), u_inp.end(), xvel_bc_data[lev].begin());
+
679  Gpu::copy(Gpu::hostToDevice, v_inp.begin(), v_inp.end(), yvel_bc_data[lev].begin());
+
680  Gpu::copy(Gpu::hostToDevice, w_inp.begin(), w_inp.end(), zvel_bc_data[lev].begin());
+
681 
+
682  // NOTE: These device vectors are passed to the PhysBC constructors when that
+
683  // class is instantiated in ERF_MakeNewArrays.cpp.
+
684  } // lev
+
685 }
Coord::z
@ z
interpolate_1d
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real interpolate_1d(const amrex::Real *alpha, const amrex::Real *beta, const amrex::Real alpha_interp, const int alpha_size)
Definition: ERF_Interpolation_1D.H:12
@@ -13963,220 +13967,220 @@

Function to define the refinement criteria based on user input

-
106 {
-
107  if (max_level > 0)
-
108  {
-
109  ParmParse pp(pp_prefix);
-
110  Vector<std::string> refinement_indicators;
-
111  pp.queryarr("refinement_indicators",refinement_indicators,0,pp.countval("refinement_indicators"));
-
112 
-
113  for (int i=0; i<refinement_indicators.size(); ++i)
-
114  {
-
115  std::string ref_prefix = pp_prefix + "." + refinement_indicators[i];
-
116 
-
117  ParmParse ppr(ref_prefix);
-
118  RealBox realbox;
-
119  int lev_for_box;
-
120 
-
121  int num_real_lo = ppr.countval("in_box_lo");
-
122  int num_indx_lo = ppr.countval("in_box_lo_indices");
+
109 {
+
110  if (max_level > 0)
+
111  {
+
112  ParmParse pp(pp_prefix);
+
113  Vector<std::string> refinement_indicators;
+
114  pp.queryarr("refinement_indicators",refinement_indicators,0,pp.countval("refinement_indicators"));
+
115 
+
116  for (int i=0; i<refinement_indicators.size(); ++i)
+
117  {
+
118  std::string ref_prefix = pp_prefix + "." + refinement_indicators[i];
+
119 
+
120  ParmParse ppr(ref_prefix);
+
121  RealBox realbox;
+
122  int lev_for_box;
123 
-
124  if ( !((num_real_lo == AMREX_SPACEDIM && num_indx_lo == 0) ||
-
125  (num_indx_lo == AMREX_SPACEDIM && num_real_lo == 0) ||
-
126  (num_indx_lo == 0 && num_real_lo == 0)) )
-
127  {
-
128  amrex::Abort("Must only specify box for refinement using real OR index space");
-
129  }
-
130 
-
131  if (num_real_lo > 0) {
-
132  std::vector<Real> rbox_lo(3), rbox_hi(3);
-
133  ppr.get("max_level",lev_for_box);
-
134  if (lev_for_box <= max_level)
-
135  {
-
136  if (n_error_buf[0] != IntVect::TheZeroVector()) {
-
137  amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
-
138  }
-
139 
-
140  const Real* plo = geom[lev_for_box].ProbLo();
-
141  const Real* phi = geom[lev_for_box].ProbHi();
+
124  int num_real_lo = ppr.countval("in_box_lo");
+
125  int num_indx_lo = ppr.countval("in_box_lo_indices");
+
126 
+
127  if ( !((num_real_lo == AMREX_SPACEDIM && num_indx_lo == 0) ||
+
128  (num_indx_lo == AMREX_SPACEDIM && num_real_lo == 0) ||
+
129  (num_indx_lo == 0 && num_real_lo == 0)) )
+
130  {
+
131  amrex::Abort("Must only specify box for refinement using real OR index space");
+
132  }
+
133 
+
134  if (num_real_lo > 0) {
+
135  std::vector<Real> rbox_lo(3), rbox_hi(3);
+
136  ppr.get("max_level",lev_for_box);
+
137  if (lev_for_box <= max_level)
+
138  {
+
139  if (n_error_buf[0] != IntVect::TheZeroVector()) {
+
140  amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
+
141  }
142 
-
143  ppr.getarr("in_box_lo",rbox_lo,0,AMREX_SPACEDIM);
-
144  ppr.getarr("in_box_hi",rbox_hi,0,AMREX_SPACEDIM);
+
143  const Real* plo = geom[lev_for_box].ProbLo();
+
144  const Real* phi = geom[lev_for_box].ProbHi();
145 
-
146  if (rbox_lo[0] < plo[0]) rbox_lo[0] = plo[0];
-
147  if (rbox_lo[1] < plo[1]) rbox_lo[1] = plo[1];
-
148  if (rbox_hi[0] > phi[0]) rbox_hi[0] = phi[0];
-
149  if (rbox_hi[1] > phi[1]) rbox_hi[1] = phi[1];
-
150 
-
151  realbox = RealBox(&(rbox_lo[0]),&(rbox_hi[0]));
-
152 
-
153  Print() << "Realbox read in and intersected laterally with domain is " << realbox << std::endl;
-
154 
-
155  num_boxes_at_level[lev_for_box] += 1;
-
156 
-
157  int ilo, jlo, klo;
-
158  int ihi, jhi, khi;
-
159  const auto* dx = geom[lev_for_box].CellSize();
-
160  ilo = static_cast<int>((rbox_lo[0] - plo[0])/dx[0]);
-
161  jlo = static_cast<int>((rbox_lo[1] - plo[1])/dx[1]);
-
162  ihi = static_cast<int>((rbox_hi[0] - plo[0])/dx[0]-1);
-
163  jhi = static_cast<int>((rbox_hi[1] - plo[1])/dx[1]-1);
-
164  if (SolverChoice::mesh_type != MeshType::ConstantDz) {
-
165  // Search for k indices corresponding to nominal grid
-
166  // AGL heights
-
167  const Box& domain = geom[lev_for_box].Domain();
-
168  klo = domain.smallEnd(2) - 1;
-
169  khi = domain.smallEnd(2) - 1;
-
170 
-
171  if (rbox_lo[2] < zlevels_stag[lev_for_box][domain.smallEnd(2)])
-
172  {
-
173  klo = domain.smallEnd(2);
-
174  }
-
175  else
-
176  {
-
177  for (int k=domain.smallEnd(2); k<=domain.bigEnd(2)+1; ++k) {
-
178  if (zlevels_stag[lev_for_box][k] > rbox_lo[2]) {
-
179  klo = k-1;
-
180  break;
-
181  }
-
182  }
-
183  }
-
184  AMREX_ASSERT(klo >= domain.smallEnd(2));
-
185 
-
186  if (rbox_hi[2] > zlevels_stag[lev_for_box][domain.bigEnd(2)+1])
-
187  {
-
188  khi = domain.bigEnd(2);
-
189  }
-
190  else
-
191  {
-
192  for (int k=klo+1; k<=domain.bigEnd(2)+1; ++k) {
-
193  if (zlevels_stag[lev_for_box][k] > rbox_hi[2]) {
-
194  khi = k-1;
-
195  break;
-
196  }
-
197  }
-
198  }
-
199  AMREX_ASSERT((khi <= domain.bigEnd(2)) && (khi > klo));
-
200 
-
201  // Need to update realbox because tagging is based on
-
202  // the initial _un_deformed grid
-
203  realbox = RealBox(plo[0]+ ilo *dx[0], plo[1]+ jlo *dx[1], plo[2]+ klo *dx[2],
-
204  plo[0]+(ihi+1)*dx[0], plo[1]+(jhi+1)*dx[1], plo[2]+(khi+1)*dx[2]);
-
205  } else {
-
206  klo = static_cast<int>((rbox_lo[2] - plo[2])/dx[2]);
-
207  khi = static_cast<int>((rbox_hi[2] - plo[2])/dx[2]-1);
-
208  }
-
209 
-
210  Box bx(IntVect(ilo,jlo,klo),IntVect(ihi,jhi,khi));
-
211  if ( (ilo%ref_ratio[lev_for_box-1][0] != 0) || ((ihi+1)%ref_ratio[lev_for_box-1][0] != 0) ||
-
212  (jlo%ref_ratio[lev_for_box-1][1] != 0) || ((jhi+1)%ref_ratio[lev_for_box-1][1] != 0) ||
-
213  (klo%ref_ratio[lev_for_box-1][2] != 0) || ((khi+1)%ref_ratio[lev_for_box-1][2] != 0) )
-
214  {
-
215  amrex::Print() << "Box : " << bx << std::endl;
-
216  amrex::Print() << "RealBox : " << realbox << std::endl;
-
217  amrex::Print() << "ilo, ihi+1, jlo, jhi+1, klo, khi+1 by ref_ratio : "
-
218  << ilo%ref_ratio[lev_for_box-1][0] << " " << (ihi+1)%ref_ratio[lev_for_box-1][0] << " "
-
219  << jlo%ref_ratio[lev_for_box-1][1] << " " << (jhi+1)%ref_ratio[lev_for_box-1][1] << " "
-
220  << klo%ref_ratio[lev_for_box-1][2] << " " << (khi+1)%ref_ratio[lev_for_box-1][2] << std::endl;
-
221  amrex::Error("Fine box is not legit with this ref_ratio");
-
222  }
-
223  boxes_at_level[lev_for_box].push_back(bx);
-
224  Print() << "Saving in 'boxes at level' as " << bx << std::endl;
-
225  } // lev
-
226  if (init_type == InitType::Real || init_type == InitType::Metgrid) {
-
227  if (num_boxes_at_level[lev_for_box] != num_files_at_level[lev_for_box]) {
-
228  amrex::Error("Number of boxes doesn't match number of input files");
-
229 
-
230  }
-
231  }
+
146  ppr.getarr("in_box_lo",rbox_lo,0,AMREX_SPACEDIM);
+
147  ppr.getarr("in_box_hi",rbox_hi,0,AMREX_SPACEDIM);
+
148 
+
149  if (rbox_lo[0] < plo[0]) rbox_lo[0] = plo[0];
+
150  if (rbox_lo[1] < plo[1]) rbox_lo[1] = plo[1];
+
151  if (rbox_hi[0] > phi[0]) rbox_hi[0] = phi[0];
+
152  if (rbox_hi[1] > phi[1]) rbox_hi[1] = phi[1];
+
153 
+
154  realbox = RealBox(&(rbox_lo[0]),&(rbox_hi[0]));
+
155 
+
156  Print() << "Realbox read in and intersected laterally with domain is " << realbox << std::endl;
+
157 
+
158  num_boxes_at_level[lev_for_box] += 1;
+
159 
+
160  int ilo, jlo, klo;
+
161  int ihi, jhi, khi;
+
162  const auto* dx = geom[lev_for_box].CellSize();
+
163  ilo = static_cast<int>((rbox_lo[0] - plo[0])/dx[0]);
+
164  jlo = static_cast<int>((rbox_lo[1] - plo[1])/dx[1]);
+
165  ihi = static_cast<int>((rbox_hi[0] - plo[0])/dx[0]-1);
+
166  jhi = static_cast<int>((rbox_hi[1] - plo[1])/dx[1]-1);
+
167  if (SolverChoice::mesh_type != MeshType::ConstantDz) {
+
168  // Search for k indices corresponding to nominal grid
+
169  // AGL heights
+
170  const Box& domain = geom[lev_for_box].Domain();
+
171  klo = domain.smallEnd(2) - 1;
+
172  khi = domain.smallEnd(2) - 1;
+
173 
+
174  if (rbox_lo[2] < zlevels_stag[lev_for_box][domain.smallEnd(2)])
+
175  {
+
176  klo = domain.smallEnd(2);
+
177  }
+
178  else
+
179  {
+
180  for (int k=domain.smallEnd(2); k<=domain.bigEnd(2)+1; ++k) {
+
181  if (zlevels_stag[lev_for_box][k] > rbox_lo[2]) {
+
182  klo = k-1;
+
183  break;
+
184  }
+
185  }
+
186  }
+
187  AMREX_ASSERT(klo >= domain.smallEnd(2));
+
188 
+
189  if (rbox_hi[2] > zlevels_stag[lev_for_box][domain.bigEnd(2)+1])
+
190  {
+
191  khi = domain.bigEnd(2);
+
192  }
+
193  else
+
194  {
+
195  for (int k=klo+1; k<=domain.bigEnd(2)+1; ++k) {
+
196  if (zlevels_stag[lev_for_box][k] > rbox_hi[2]) {
+
197  khi = k-1;
+
198  break;
+
199  }
+
200  }
+
201  }
+
202  AMREX_ASSERT((khi <= domain.bigEnd(2)) && (khi > klo));
+
203 
+
204  // Need to update realbox because tagging is based on
+
205  // the initial _un_deformed grid
+
206  realbox = RealBox(plo[0]+ ilo *dx[0], plo[1]+ jlo *dx[1], plo[2]+ klo *dx[2],
+
207  plo[0]+(ihi+1)*dx[0], plo[1]+(jhi+1)*dx[1], plo[2]+(khi+1)*dx[2]);
+
208  } else {
+
209  klo = static_cast<int>((rbox_lo[2] - plo[2])/dx[2]);
+
210  khi = static_cast<int>((rbox_hi[2] - plo[2])/dx[2]-1);
+
211  }
+
212 
+
213  Box bx(IntVect(ilo,jlo,klo),IntVect(ihi,jhi,khi));
+
214  if ( (ilo%ref_ratio[lev_for_box-1][0] != 0) || ((ihi+1)%ref_ratio[lev_for_box-1][0] != 0) ||
+
215  (jlo%ref_ratio[lev_for_box-1][1] != 0) || ((jhi+1)%ref_ratio[lev_for_box-1][1] != 0) ||
+
216  (klo%ref_ratio[lev_for_box-1][2] != 0) || ((khi+1)%ref_ratio[lev_for_box-1][2] != 0) )
+
217  {
+
218  amrex::Print() << "Box : " << bx << std::endl;
+
219  amrex::Print() << "RealBox : " << realbox << std::endl;
+
220  amrex::Print() << "ilo, ihi+1, jlo, jhi+1, klo, khi+1 by ref_ratio : "
+
221  << ilo%ref_ratio[lev_for_box-1][0] << " " << (ihi+1)%ref_ratio[lev_for_box-1][0] << " "
+
222  << jlo%ref_ratio[lev_for_box-1][1] << " " << (jhi+1)%ref_ratio[lev_for_box-1][1] << " "
+
223  << klo%ref_ratio[lev_for_box-1][2] << " " << (khi+1)%ref_ratio[lev_for_box-1][2] << std::endl;
+
224  amrex::Error("Fine box is not legit with this ref_ratio");
+
225  }
+
226  boxes_at_level[lev_for_box].push_back(bx);
+
227  Print() << "Saving in 'boxes at level' as " << bx << std::endl;
+
228  } // lev
+
229  if (init_type == InitType::Real || init_type == InitType::Metgrid) {
+
230  if (num_boxes_at_level[lev_for_box] != num_files_at_level[lev_for_box]) {
+
231  amrex::Error("Number of boxes doesn't match number of input files");
232 
-
233  } else if (num_indx_lo > 0) {
-
234 
-
235  std::vector<int> box_lo(3), box_hi(3);
-
236  ppr.get("max_level",lev_for_box);
-
237  if (lev_for_box <= max_level)
-
238  {
-
239  if (n_error_buf[0] != IntVect::TheZeroVector()) {
-
240  amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
-
241  }
-
242 
-
243  ppr.getarr("in_box_lo_indices",box_lo,0,AMREX_SPACEDIM);
-
244  ppr.getarr("in_box_hi_indices",box_hi,0,AMREX_SPACEDIM);
+
233  }
+
234  }
+
235 
+
236  } else if (num_indx_lo > 0) {
+
237 
+
238  std::vector<int> box_lo(3), box_hi(3);
+
239  ppr.get("max_level",lev_for_box);
+
240  if (lev_for_box <= max_level)
+
241  {
+
242  if (n_error_buf[0] != IntVect::TheZeroVector()) {
+
243  amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
+
244  }
245 
-
246  Box bx(IntVect(box_lo[0],box_lo[1],box_lo[2]),IntVect(box_hi[0],box_hi[1],box_hi[2]));
-
247  amrex::Print() << "BOX " << bx << std::endl;
+
246  ppr.getarr("in_box_lo_indices",box_lo,0,AMREX_SPACEDIM);
+
247  ppr.getarr("in_box_hi_indices",box_hi,0,AMREX_SPACEDIM);
248 
-
249  const auto* dx = geom[lev_for_box].CellSize();
-
250  const Real* plo = geom[lev_for_box].ProbLo();
-
251  realbox = RealBox(plo[0]+ box_lo[0] *dx[0], plo[1]+ box_lo[1] *dx[1], plo[2]+ box_lo[2] *dx[2],
-
252  plo[0]+(box_hi[0]+1)*dx[0], plo[1]+(box_hi[1]+1)*dx[1], plo[2]+(box_hi[2]+1)*dx[2]);
-
253 
-
254  Print() << "Reading " << bx << " at level " << lev_for_box << std::endl;
-
255  num_boxes_at_level[lev_for_box] += 1;
+
249  Box bx(IntVect(box_lo[0],box_lo[1],box_lo[2]),IntVect(box_hi[0],box_hi[1],box_hi[2]));
+
250  amrex::Print() << "BOX " << bx << std::endl;
+
251 
+
252  const auto* dx = geom[lev_for_box].CellSize();
+
253  const Real* plo = geom[lev_for_box].ProbLo();
+
254  realbox = RealBox(plo[0]+ box_lo[0] *dx[0], plo[1]+ box_lo[1] *dx[1], plo[2]+ box_lo[2] *dx[2],
+
255  plo[0]+(box_hi[0]+1)*dx[0], plo[1]+(box_hi[1]+1)*dx[1], plo[2]+(box_hi[2]+1)*dx[2]);
256 
-
257  if ( (box_lo[0]%ref_ratio[lev_for_box-1][0] != 0) || ((box_hi[0]+1)%ref_ratio[lev_for_box-1][0] != 0) ||
-
258  (box_lo[1]%ref_ratio[lev_for_box-1][1] != 0) || ((box_hi[1]+1)%ref_ratio[lev_for_box-1][1] != 0) ||
-
259  (box_lo[2]%ref_ratio[lev_for_box-1][2] != 0) || ((box_hi[2]+1)%ref_ratio[lev_for_box-1][2] != 0) )
-
260  amrex::Error("Fine box is not legit with this ref_ratio");
-
261  boxes_at_level[lev_for_box].push_back(bx);
-
262  Print() << "Saving in 'boxes at level' as " << bx << std::endl;
-
263  } // lev
-
264  if (init_type == InitType::Real || init_type == InitType::Metgrid) {
-
265  if (num_boxes_at_level[lev_for_box] != num_files_at_level[lev_for_box]) {
-
266  amrex::Error("Number of boxes doesn't match number of input files");
-
267 
-
268  }
-
269  }
-
270  }
-
271 
-
272  AMRErrorTagInfo info;
-
273 
-
274  if (realbox.ok()) {
-
275  info.SetRealBox(realbox);
-
276  }
-
277  if (ppr.countval("start_time") > 0) {
-
278  Real ref_min_time; ppr.get("start_time",ref_min_time);
-
279  info.SetMinTime(ref_min_time);
-
280  }
-
281  if (ppr.countval("end_time") > 0) {
-
282  Real ref_max_time; ppr.get("end_time",ref_max_time);
-
283  info.SetMaxTime(ref_max_time);
-
284  }
-
285  if (ppr.countval("max_level") > 0) {
-
286  int ref_max_level; ppr.get("max_level",ref_max_level);
-
287  info.SetMaxLevel(ref_max_level);
-
288  }
-
289 
-
290  if (ppr.countval("value_greater")) {
-
291  int num_val = ppr.countval("value_greater");
-
292  Vector<Real> value(num_val);
-
293  ppr.getarr("value_greater",value,0,num_val);
-
294  std::string field; ppr.get("field_name",field);
-
295  ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::GREATER,field,info));
-
296  }
-
297  else if (ppr.countval("value_less")) {
-
298  int num_val = ppr.countval("value_less");
-
299  Vector<Real> value(num_val);
-
300  ppr.getarr("value_less",value,0,num_val);
-
301  std::string field; ppr.get("field_name",field);
-
302  ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::LESS,field,info));
-
303  }
-
304  else if (ppr.countval("adjacent_difference_greater")) {
-
305  int num_val = ppr.countval("adjacent_difference_greater");
-
306  Vector<Real> value(num_val);
-
307  ppr.getarr("adjacent_difference_greater",value,0,num_val);
-
308  std::string field; ppr.get("field_name",field);
-
309  ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::GRAD,field,info));
-
310  }
-
311  else if (realbox.ok())
-
312  {
-
313  ref_tags.push_back(AMRErrorTag(info));
-
314  } else {
-
315  Abort(std::string("Unrecognized refinement indicator for " + refinement_indicators[i]).c_str());
-
316  }
-
317  } // loop over criteria
-
318  } // if max_level > 0
-
319 }
+
257  Print() << "Reading " << bx << " at level " << lev_for_box << std::endl;
+
258  num_boxes_at_level[lev_for_box] += 1;
+
259 
+
260  if ( (box_lo[0]%ref_ratio[lev_for_box-1][0] != 0) || ((box_hi[0]+1)%ref_ratio[lev_for_box-1][0] != 0) ||
+
261  (box_lo[1]%ref_ratio[lev_for_box-1][1] != 0) || ((box_hi[1]+1)%ref_ratio[lev_for_box-1][1] != 0) ||
+
262  (box_lo[2]%ref_ratio[lev_for_box-1][2] != 0) || ((box_hi[2]+1)%ref_ratio[lev_for_box-1][2] != 0) )
+
263  amrex::Error("Fine box is not legit with this ref_ratio");
+
264  boxes_at_level[lev_for_box].push_back(bx);
+
265  Print() << "Saving in 'boxes at level' as " << bx << std::endl;
+
266  } // lev
+
267  if (init_type == InitType::Real || init_type == InitType::Metgrid) {
+
268  if (num_boxes_at_level[lev_for_box] != num_files_at_level[lev_for_box]) {
+
269  amrex::Error("Number of boxes doesn't match number of input files");
+
270 
+
271  }
+
272  }
+
273  }
+
274 
+
275  AMRErrorTagInfo info;
+
276 
+
277  if (realbox.ok()) {
+
278  info.SetRealBox(realbox);
+
279  }
+
280  if (ppr.countval("start_time") > 0) {
+
281  Real ref_min_time; ppr.get("start_time",ref_min_time);
+
282  info.SetMinTime(ref_min_time);
+
283  }
+
284  if (ppr.countval("end_time") > 0) {
+
285  Real ref_max_time; ppr.get("end_time",ref_max_time);
+
286  info.SetMaxTime(ref_max_time);
+
287  }
+
288  if (ppr.countval("max_level") > 0) {
+
289  int ref_max_level; ppr.get("max_level",ref_max_level);
+
290  info.SetMaxLevel(ref_max_level);
+
291  }
+
292 
+
293  if (ppr.countval("value_greater")) {
+
294  int num_val = ppr.countval("value_greater");
+
295  Vector<Real> value(num_val);
+
296  ppr.getarr("value_greater",value,0,num_val);
+
297  std::string field; ppr.get("field_name",field);
+
298  ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::GREATER,field,info));
+
299  }
+
300  else if (ppr.countval("value_less")) {
+
301  int num_val = ppr.countval("value_less");
+
302  Vector<Real> value(num_val);
+
303  ppr.getarr("value_less",value,0,num_val);
+
304  std::string field; ppr.get("field_name",field);
+
305  ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::LESS,field,info));
+
306  }
+
307  else if (ppr.countval("adjacent_difference_greater")) {
+
308  int num_val = ppr.countval("adjacent_difference_greater");
+
309  Vector<Real> value(num_val);
+
310  ppr.getarr("adjacent_difference_greater",value,0,num_val);
+
311  std::string field; ppr.get("field_name",field);
+
312  ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::GRAD,field,info));
+
313  }
+
314  else if (realbox.ok())
+
315  {
+
316  ref_tags.push_back(AMRErrorTag(info));
+
317  } else {
+
318  Abort(std::string("Unrecognized refinement indicator for " + refinement_indicators[i]).c_str());
+
319  }
+
320  } // loop over criteria
+
321  } // if max_level > 0
+
322 }
Here is the call graph for this function: