ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ReadBndryPlanes Class Reference

#include <ERF_ReadBndryPlanes.H>

Collaboration diagram for ReadBndryPlanes:

Public Member Functions

 ReadBndryPlanes (const amrex::Geometry &geom, const amrex::Real &rdOcp_in)
 
void define_level_data (int lev)
 
void read_time_file ()
 
void read_input_files (amrex::Real time, amrex::Real dt, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_extdir_vals)
 
void read_file (int idx, amrex::Vector< std::unique_ptr< PlaneVector >> &data_to_fill, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_extdir_vals)
 
amrex::Vector< std::unique_ptr< PlaneVector > > & interp_in_time (const amrex::Real &time)
 
amrex::Vector< std::unique_ptr< PlaneVector > > & get_tendency (const amrex::Real &time)
 
amrex::Real tinterp () const
 
int ingested_velocity () const
 
int ingested_theta () const
 
int ingested_density () const
 
int ingested_scalar () const
 
int ingested_q1 () const
 
int ingested_q2 () const
 
int ingested_KE () const
 

Private Attributes

amrex::Real m_tn
 The times for which we currently have data. More...
 
amrex::Real m_tnp1
 
amrex::Real m_tnp2
 
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_n
 Data at time m_tn. More...
 
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_np1
 Data at time m_tnp1. More...
 
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_np2
 Data at time m_tnp2. More...
 
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_interp
 Data interpolated to the time requested. More...
 
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_tendency
 Tendency between the n and np1 data. More...
 
amrex::Real m_tinterp {-one}
 Time for plane at interpolation. More...
 
amrex::Geometry m_geom
 Geometry at level 0. More...
 
std::string m_filename {""}
 File name for IO. More...
 
std::string m_time_file {""}
 File name for file holding timesteps and times. More...
 
amrex::Vector< amrex::Realm_in_times
 The timesteps / times that we read from time.dat. More...
 
amrex::Vector< int > m_in_timesteps
 
amrex::Vector< std::string > m_var_names
 Variables to be read in. More...
 
int m_in_rad = 1
 Controls extents on native bndry output. More...
 
const int m_out_rad = 1
 
const int m_extent_rad = 0
 
bool m_use_real_bcs = false
 Are real BCs being used? More...
 
const amrex::Real m_rdOcp
 R_d/c_p is needed for reading boundary files. More...
 
int is_velocity_read
 
int is_density_read
 
int is_temperature_read
 
int is_theta_read
 
int is_scalar_read
 
int is_q1_read
 
int is_q2_read
 
int is_KE_read
 
int last_file_read
 

Detailed Description

Collection of data structures and operations for reading data

This class contains the inlet data structures and operations to read and interpolate inflow data.

Constructor & Destructor Documentation

◆ ReadBndryPlanes()

ReadBndryPlanes::ReadBndryPlanes ( const amrex::Geometry &  geom,
const amrex::Real rdOcp_in 
)
explicit

ReadBndryPlanes class constructor. Handles initialization from inputs file parameters.

Parameters
geomGeometry for the domain
rdOcp_inReal constant for the Rhydberg constant ($R_d$) divided by the specific heat at constant pressure ($c_p$)
202 :
203  m_geom(geom),
204  m_rdOcp(rdOcp_in)
205 {
206  ParmParse pp("erf");
207 
208  // Get the radius inside the domain
209  pp.query("in_rad",m_in_rad);
210 
211  // Are we using real bcs?
212  pp.query("use_real_bcs", m_use_real_bcs);
213 
214  last_file_read = -1;
215 
216  m_tinterp = -one;
217 
218  // What folder will the time series of planes be read from
219  pp.get("bndry_file", m_filename);
220 
221  is_velocity_read = 0;
222  is_density_read = 0;
224  is_theta_read = 0;
225  is_scalar_read = 0;
226  is_q1_read = 0;
227  is_q2_read = 0;
228  is_KE_read = 0;
229 
230  if (pp.contains("bndry_input_var_names"))
231  {
232  int num_vars = pp.countval("bndry_input_var_names");
233  m_var_names.resize(num_vars);
234  pp.queryarr("bndry_input_var_names",m_var_names,0,num_vars);
235  for (int i = 0; i < m_var_names.size(); i++) {
236  if (m_var_names[i] == "velocity") is_velocity_read = 1;
237  if (m_var_names[i] == "density") is_density_read = 1;
238  if (m_var_names[i] == "temperature") is_temperature_read = 1;
239  if (m_var_names[i] == "theta") is_theta_read = 1;
240  if (m_var_names[i] == "scalar") is_scalar_read = 1;
241  if (m_var_names[i] == "qv") is_q1_read = 1;
242  if (m_var_names[i] == "qc") is_q2_read = 1;
243  if (m_var_names[i] == "ke") is_KE_read = 1;
244  }
245  }
246 
247  // time.dat will be in the same folder as the time series of data
248  m_time_file = m_filename + "/time.dat";
249 
250  // each pointer (at at given time) has 6 components, one for each orientation
251  // TODO: we really only need 4 not 6
252  int size = 2*AMREX_SPACEDIM;
253  m_data_n.resize(size);
254  m_data_np1.resize(size);
255  m_data_np2.resize(size);
256  m_data_interp.resize(size);
257  m_data_tendency.resize(size);
258 }
constexpr amrex::Real one
Definition: ERF_Constants.H:9
ParmParse pp("prob")
int is_velocity_read
Definition: ERF_ReadBndryPlanes.H:104
int is_q2_read
Definition: ERF_ReadBndryPlanes.H:110
int is_theta_read
Definition: ERF_ReadBndryPlanes.H:107
bool m_use_real_bcs
Are real BCs being used?
Definition: ERF_ReadBndryPlanes.H:99
std::string m_filename
File name for IO.
Definition: ERF_ReadBndryPlanes.H:81
amrex::Real m_tinterp
Time for plane at interpolation.
Definition: ERF_ReadBndryPlanes.H:75
int is_temperature_read
Definition: ERF_ReadBndryPlanes.H:106
int is_density_read
Definition: ERF_ReadBndryPlanes.H:105
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_np2
Data at time m_tnp2.
Definition: ERF_ReadBndryPlanes.H:66
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_interp
Data interpolated to the time requested.
Definition: ERF_ReadBndryPlanes.H:69
int last_file_read
Definition: ERF_ReadBndryPlanes.H:113
int is_KE_read
Definition: ERF_ReadBndryPlanes.H:111
const amrex::Real m_rdOcp
R_d/c_p is needed for reading boundary files.
Definition: ERF_ReadBndryPlanes.H:102
std::string m_time_file
File name for file holding timesteps and times.
Definition: ERF_ReadBndryPlanes.H:84
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_tendency
Tendency between the n and np1 data.
Definition: ERF_ReadBndryPlanes.H:72
amrex::Vector< std::string > m_var_names
Variables to be read in.
Definition: ERF_ReadBndryPlanes.H:91
int is_scalar_read
Definition: ERF_ReadBndryPlanes.H:108
int is_q1_read
Definition: ERF_ReadBndryPlanes.H:109
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_n
Data at time m_tn.
Definition: ERF_ReadBndryPlanes.H:60
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_np1
Data at time m_tnp1.
Definition: ERF_ReadBndryPlanes.H:63
int m_in_rad
Controls extents on native bndry output.
Definition: ERF_ReadBndryPlanes.H:94
amrex::Geometry m_geom
Geometry at level 0.
Definition: ERF_ReadBndryPlanes.H:78
Here is the call graph for this function:

Member Function Documentation

◆ define_level_data()

void ReadBndryPlanes::define_level_data ( int  lev)

Function in ReadBndryPlanes class for allocating space for the boundary plane data ERF will need.

44 {
45  Print() << "ReadBndryPlanes::define_level_data" << std::endl;
46  // *********************************************************
47  // Allocate space for all of the boundary planes we may need
48  // *********************************************************
49  int ncomp = BCVars::NumTypes;
50  const Box& domain = m_geom.Domain();
51  for (OrientationIter oit; oit != nullptr; ++oit) {
52  auto ori = oit();
53  if (ori.coordDir() < 2) {
54 
55  m_data_n[ori] = std::make_unique<PlaneVector>();
56  m_data_np1[ori] = std::make_unique<PlaneVector>();
57  m_data_np2[ori] = std::make_unique<PlaneVector>();
58  m_data_interp[ori] = std::make_unique<PlaneVector>();
59  m_data_tendency[ori] = std::make_unique<PlaneVector>();
60 
61  const auto& lo = domain.loVect();
62  const auto& hi = domain.hiVect();
63 
64  IntVect plo(lo);
65  IntVect phi(hi);
66  const int normal = ori.coordDir();
67  plo[normal] = ori.isHigh() ? hi[normal] - (m_in_rad - 1) : -m_out_rad;
68  phi[normal] = ori.isHigh() ? hi[normal] + (m_out_rad ) : (m_in_rad - 1);
69  const Box pbx(plo, phi);
70  m_data_n[ori]->push_back(FArrayBox(pbx, ncomp));
71  m_data_np1[ori]->push_back(FArrayBox(pbx, ncomp));
72  m_data_np2[ori]->push_back(FArrayBox(pbx, ncomp));
73  m_data_interp[ori]->push_back(FArrayBox(pbx, ncomp));
74  m_data_tendency[ori]->push_back(FArrayBox(pbx, ncomp));
75  }
76  }
77 }
const int m_out_rad
Definition: ERF_ReadBndryPlanes.H:95
@ NumTypes
Definition: ERF_IndexDefines.H:105

Referenced by read_time_file().

Here is the caller graph for this function:

◆ get_tendency()

Vector< std::unique_ptr< PlaneVector > > & ReadBndryPlanes::get_tendency ( const amrex::Real time)

Function in ReadBndryPlanes class for interpolating boundary data in time.

Parameters
timeConstant specifying the time for interpolation
145 {
146  AMREX_ALWAYS_ASSERT(m_tn <= time && time <= m_tnp2);
147 
148  if (time < m_tnp1) {
149  Real idt = Real(1.0) / (m_tnp1 - m_tn);
150  for (OrientationIter oit; oit != nullptr; ++oit) {
151  auto ori = oit();
152  if (ori.coordDir() < 2) {
153  const int nlevels = static_cast<int>(m_data_n[ori]->size());
154  for (int lev = 0; lev < nlevels; ++lev) {
155  auto& fabt = (*m_data_tendency[ori])[lev];
156  Box bx = fabt.box();
157  int ncomp = fabt.nComp();
158 
159  const auto& datt = fabt.array();
160  const auto& datn = (*m_data_n[ori])[lev].array();
161  const auto& datnp1 = (*m_data_np1[ori])[lev].array();
162  ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
163  {
164  datt(i,j,k,n) = (datnp1(i,j,k,n) - datn(i,j,k,n)) * idt;
165  });
166  }
167  }
168  }
169  } else {
170  Real idt = Real(1.0) / (m_tnp2 - m_tnp1);
171  for (OrientationIter oit; oit != nullptr; ++oit) {
172  auto ori = oit();
173  if (ori.coordDir() < 2) {
174  const int nlevels = static_cast<int>(m_data_n[ori]->size());
175  for (int lev = 0; lev < nlevels; ++lev) {
176  auto& fabt = (*m_data_tendency[ori])[lev];
177  Box bx = fabt.box();
178  int ncomp = fabt.nComp();
179 
180  const auto& datt = fabt.array();
181  const auto& datnp1 = (*m_data_np1[ori])[lev].array();
182  const auto& datnp2 = (*m_data_np2[ori])[lev].array();
183  ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
184  {
185  datt(i,j,k,n) = (datnp2(i,j,k,n) - datnp1(i,j,k,n)) * idt;
186  });
187  }
188  }
189  }
190  }
191 
192  return m_data_tendency;
193 }
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor(grown_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
amrex::Real Real
Definition: ERF_ShocInterface.H:19
amrex::Real m_tnp1
Definition: ERF_ReadBndryPlanes.H:56
amrex::Real m_tnp2
Definition: ERF_ReadBndryPlanes.H:57
amrex::Real m_tn
The times for which we currently have data.
Definition: ERF_ReadBndryPlanes.H:55
Here is the call graph for this function:

◆ ingested_density()

int ReadBndryPlanes::ingested_density ( ) const
inline
46 {return is_density_read;}

Referenced by read_file().

Here is the caller graph for this function:

◆ ingested_KE()

int ReadBndryPlanes::ingested_KE ( ) const
inline
50 {return is_KE_read;}

◆ ingested_q1()

int ReadBndryPlanes::ingested_q1 ( ) const
inline
48 {return is_q1_read;}

◆ ingested_q2()

int ReadBndryPlanes::ingested_q2 ( ) const
inline
49 {return is_q2_read;}

◆ ingested_scalar()

int ReadBndryPlanes::ingested_scalar ( ) const
inline
47 {return is_scalar_read;}

◆ ingested_theta()

int ReadBndryPlanes::ingested_theta ( ) const
inline

◆ ingested_velocity()

int ReadBndryPlanes::ingested_velocity ( ) const
inline
44 {return is_velocity_read;}

◆ interp_in_time()

Vector< std::unique_ptr< PlaneVector > > & ReadBndryPlanes::interp_in_time ( const amrex::Real time)

Function in ReadBndryPlanes class for interpolating boundary data in time.

Parameters
timeConstant specifying the time for interpolation
87 {
88  AMREX_ALWAYS_ASSERT(m_tn <= time && time <= m_tnp2);
89 
90  //Print() << "interp_in_time at time " << time << " given " << m_tn << " " << m_tnp1 << " " << m_tnp2 << std::endl;
91  //Print() << "m_tinterp " << m_tinterp << std::endl;
92 
93  if (time == m_tinterp) {
94  // We have already interpolated to this time
95  return m_data_interp;
96 
97  } else {
98 
99  // We must now interpolate to a new time
100  m_tinterp = time;
101 
102  if (time < m_tnp1) {
103  for (OrientationIter oit; oit != nullptr; ++oit) {
104  auto ori = oit();
105  if (ori.coordDir() < 2) {
106  const int nlevels = static_cast<int>(m_data_n[ori]->size());
107  for (int lev = 0; lev < nlevels; ++lev) {
108  const auto& datn = (*m_data_n[ori])[lev];
109  const auto& datnp1 = (*m_data_np1[ori])[lev];
110  auto& dati = (*m_data_interp[ori])[lev];
111  dati.linInterp<RunOn::Device>(datn, 0, datnp1, 0,
113  datn.box(), 0, dati.nComp());
114  }
115  }
116  }
117  } else {
118  for (OrientationIter oit; oit != nullptr; ++oit) {
119  auto ori = oit();
120  if (ori.coordDir() < 2) {
121  const int nlevels = static_cast<int>(m_data_n[ori]->size());
122  for (int lev = 0; lev < nlevels; ++lev) {
123  const auto& datnp1 = (*m_data_np1[ori])[lev];
124  const auto& datnp2 = (*m_data_np2[ori])[lev];
125  auto& dati = (*m_data_interp[ori])[lev];
126  dati.linInterp<RunOn::Device>(datnp1, 0, datnp2, 0,
128  datnp1.box(), 0, dati.nComp());
129  }
130  }
131  }
132  }
133  }
134  return m_data_interp;
135 }
Here is the call graph for this function:

◆ read_file()

void ReadBndryPlanes::read_file ( int  idx,
amrex::Vector< std::unique_ptr< PlaneVector >> &  data_to_fill,
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max m_bc_extdir_vals 
)

Function in ReadBndryPlanes to read boundary data for each face and variable from files.

Parameters
idxSpecifies the index corresponding to the timestep we want
data_to_fillContainer for face data on boundaries
m_bc_extdir_valsContainer storing the external dirichlet boundary conditions we are reading from the input files
411 {
412  if (idx >= m_in_timesteps.size()) {
413  Print() << "Asking for index " << idx << " but m_in_timesteps only has size " << m_in_timesteps.size() << std::endl;
414  Abort();
415  }
416  const int t_step = m_in_timesteps[idx];
417  const std::string chkname1 = m_filename + Concatenate("/bndry_output", t_step);
418 
419  const std::string level_prefix = "Level_";
420  const int lev = 0;
421 
422  const Box& domain = m_geom.Domain();
423  BoxArray ba(domain);
424  DistributionMapping dm{ba};
425 
426  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NBCVAR_max> l_bc_extdir_vals_d;
427 
428  for (int i = 0; i < BCVars::NumTypes; i++)
429  {
430  for (OrientationIter oit; oit != nullptr; ++oit) {
431  auto ori = oit();
432  l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[i][ori];
433  }
434  }
435 
436  int n_for_density = -1;
437  for (int i = 0; i < m_var_names.size(); i++)
438  {
439  if (m_var_names[i] == "density") n_for_density = i;
440  }
441 
442  // We need to initialize all the components because we may not fill all of them from files,
443  // but the loop in the interpolate routine goes over all the components anyway
444  int ncomp_for_bc = BCVars::NumTypes;
445  for (OrientationIter oit; oit != nullptr; ++oit) {
446  auto ori = oit();
447  if (ori.coordDir() < 2) {
448  FArrayBox& d = (*data_to_fill[ori])[lev];
449  const auto& bx = d.box();
450  Array4<Real> d_arr = d.array();
451  ParallelFor(
452  bx, ncomp_for_bc, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
453  d_arr(i,j,k,n) = zero;
454  });
455  }
456  }
457 
458  // Read density for primitive to conserved conversions
459  std::string filenamer = MultiFabFileFullPrefix(lev, chkname1, level_prefix, "density");
460  BndryRegister bndry_r(ba, dm, m_in_rad, m_out_rad, m_extent_rad, 1);
461  bndry_r.setVal(bogus_large_value);
462  for (OrientationIter oit; oit != nullptr; ++oit) {
463  auto ori = oit();
464  if (ori.coordDir() < 2) {
465  std::string facenamer = Concatenate(filenamer + '_', ori, 1);
466  bndry_r[ori].read(facenamer);
467  }
468  }
469 
470  // Expose for GPU
471  bool real_bcs = m_use_real_bcs;
472 
473  for (int ivar = 0; ivar < m_var_names.size(); ivar++)
474  {
475  std::string var_name = m_var_names[ivar];
476 
477  std::string filename1 = MultiFabFileFullPrefix(lev, chkname1, level_prefix, var_name);
478 
479  int ncomp;
480  if (var_name == "velocity") {
481  ncomp = AMREX_SPACEDIM;
482  } else {
483  ncomp = 1;
484  }
485 
486  int n_offset;
487  if (var_name == "density") n_offset = BCVars::Rho_bc_comp;
488  if (var_name == "theta") n_offset = BCVars::RhoTheta_bc_comp;
489  if (var_name == "temperature") n_offset = BCVars::RhoTheta_bc_comp;
490  if (var_name == "ke") n_offset = BCVars::RhoKE_bc_comp;
491  if (var_name == "scalar") n_offset = BCVars::RhoScalar_bc_comp;
492  if (var_name == "qv") n_offset = BCVars::RhoQ1_bc_comp;
493  if (var_name == "qc") n_offset = BCVars::RhoQ2_bc_comp;
494  if (var_name == "velocity") n_offset = BCVars::xvel_bc;
495 
496  // Print() << "Reading " << chkname1 << " for variable " << var_name << " with n_offset == " << n_offset << std::endl;
497 
498  BndryRegister bndry(ba, dm, m_in_rad, m_out_rad, m_extent_rad, ncomp);
499  bndry.setVal(bogus_large_value);
500 
501  // *********************************************************
502  // Read in the BndryReg for all non-z faces
503  // *********************************************************
504  for (OrientationIter oit; oit != nullptr; ++oit) {
505  auto ori = oit();
506  if (ori.coordDir() < 2) {
507 
508  std::string facename1 = Concatenate(filename1 + '_', ori, 1);
509  bndry[ori].read(facename1);
510 
511  int normal = ori.coordDir();
512  IntVect v_offset = offset(ori.faceDir(), normal);
513  if (real_bcs) { v_offset = IntVect(0); }
514 
515  const auto& bbx = (*data_to_fill[ori])[lev].box();
516 
517  // *********************************************************
518  // Copy from the BndryReg into a MultiFab then use copyTo
519  // to write from the MultiFab to a single FAB for each face
520  // *********************************************************
521  MultiFab bndryMF(
522  bndry[ori].boxArray(), bndry[ori].DistributionMap(),
523  ncomp, 0, MFInfo());
524 
525  for (MFIter mfi(bndryMF); mfi.isValid(); ++mfi) {
526 
527  const auto& vbx = mfi.validbox();
528  const auto& bndry_read_arr = bndry[ori].array(mfi);
529  const auto& bndry_read_r_arr = bndry_r[ori].array(mfi);
530  const auto& bndry_mf_arr = bndryMF.array(mfi);
531 
532  const auto& bx = bbx & vbx;
533  if (bx.isEmpty()) {
534  continue;
535  }
536 
537  // Split the 2-cell-thick working box into ghost and interior
538  // slots so the (i+v_offset) neighbor access stays in-bounds.
539  // Both slots are filled with the same face-averaged Dirichlet
540  // value; the interior slot just needs the opposite neighbor.
541  Box bx_ghost = bx;
542  Box bx_int = bx;
543  if (ori.isLow()) {
544  bx_ghost.setBig (normal, domain.smallEnd(normal) - 1);
545  bx_int .setSmall(normal, domain.smallEnd(normal));
546  } else {
547  bx_ghost.setSmall(normal, domain.bigEnd(normal) + 1);
548  bx_int .setBig (normal, domain.bigEnd(normal));
549  }
550  const IntVect v_offset_int = -v_offset;
551 
552  // We average the two cell-centered data points in the normal direction
553  // to define a Dirichlet value on the face itself.
554 
555  // This is the scalars -- they all get multiplied by rho, and in the case of
556  // reading in temperature, we must convert to theta first
557  Real rdOcp = m_rdOcp;
558  if (n_for_density >= 0) {
559  if (var_name == "temperature") {
560  ParallelFor(
561  bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
562  Real R1 = bndry_read_r_arr(i, j, k, 0);
563  Real R2 = bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
564  Real T1 = bndry_read_arr(i, j, k, 0);
565  Real T2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
566  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
567  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
568  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
569  myhalf * (R1*Th1 + R2*Th2);
570  });
571  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
572  var_name == "qv" || var_name == "qc") {
573  ParallelFor(
574  bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
575  Real R1 = bndry_read_r_arr(i, j, k, 0);
576  Real R2 = bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
577  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
578  myhalf * ( R1 * bndry_read_arr(i, j, k, 0) +
579  R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
580  });
581  } else if (var_name == "density") {
582  ParallelFor(
583  bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
584  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
585  myhalf * ( bndry_read_arr(i, j, k, 0) +
586  bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
587  });
588  }
589  } else if (!ingested_density()) {
590  if (var_name == "temperature") {
591  ParallelFor(
592  bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
593  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
594  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
595  Real T1 = bndry_read_arr(i, j, k, 0);
596  Real T2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0);
597  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
598  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
599  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
600  myhalf * (R1*Th1 + R2*Th2);
601  });
602  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
603  var_name == "qv" || var_name == "qc") {
604  ParallelFor(
605  bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
606  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
607  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
608  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
609  myhalf * (R1 * bndry_read_arr(i, j, k, 0) +
610  R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
611  });
612  }
613  }
614 
615  // This is velocity
616  if (var_name == "velocity") {
617  ParallelFor(
618  bx_ghost, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
619  bndry_mf_arr(i, j, k, n) = (real_bcs) ? bndry_read_arr(i, j, k, n) :
620  myhalf * (bndry_read_arr(i, j, k, n) +
621  bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], n));
622  });
623  }
624 
625  // --- interior-slot fill (same face-averaged Dirichlet value;
626  // neighbor offset is flipped because the "other cell"
627  // is now on the opposite side of the boundary face) ---
628  if (n_for_density >= 0) {
629  if (var_name == "temperature") {
630  ParallelFor(
631  bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
632  Real R1 = bndry_read_r_arr(i, j, k, 0);
633  Real R2 = bndry_read_r_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
634  Real T1 = bndry_read_arr(i, j, k, 0);
635  Real T2 = bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
636  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
637  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
638  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
639  myhalf * (R1*Th1 + R2*Th2);
640  });
641  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
642  var_name == "qv" || var_name == "qc") {
643  ParallelFor(
644  bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
645  Real R1 = bndry_read_r_arr(i, j, k, 0);
646  Real R2 = bndry_read_r_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
647  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
648  myhalf * ( R1 * bndry_read_arr(i, j, k, 0) +
649  R2 * bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0));
650  });
651  } else if (var_name == "density") {
652  ParallelFor(
653  bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
654  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
655  myhalf * ( bndry_read_arr(i, j, k, 0) +
656  bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0));
657  });
658  }
659  } else if (!ingested_density()) {
660  if (var_name == "temperature") {
661  ParallelFor(
662  bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
663  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
664  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
665  Real T1 = bndry_read_arr(i, j, k, 0);
666  Real T2 = bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0);
667  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
668  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
669  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
670  myhalf * (R1*Th1 + R2*Th2);
671  });
672  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
673  var_name == "qv" || var_name == "qc") {
674  ParallelFor(
675  bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
676  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
677  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
678  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
679  myhalf * (R1 * bndry_read_arr(i, j, k, 0) +
680  R2 * bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0));
681  });
682  }
683  }
684 
685  if (var_name == "velocity") {
686  ParallelFor(
687  bx_int, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
688  bndry_mf_arr(i, j, k, n) = (real_bcs) ? bndry_read_arr(i, j, k, n) :
689  myhalf * (bndry_read_arr(i, j, k, n) +
690  bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], n));
691  });
692  }
693 
694  } // mfi
695  bndryMF.copyTo((*data_to_fill[ori])[lev], 0, n_offset, ncomp);
696  } // coordDir < 2
697  } // ori
698  } // var_name
699 }
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:26
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:64
#define NBCVAR_max
Definition: ERF_IndexDefines.H:29
const Real rdOcp
Definition: ERF_InitCustomPert_Bomex.H:16
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int idx(int i, int j, int k, int nx, int ny)
Definition: ERF_InitForEnsemble.cpp:287
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
amrex::Vector< int > m_in_timesteps
Definition: ERF_ReadBndryPlanes.H:88
int ingested_density() const
Definition: ERF_ReadBndryPlanes.H:46
const int m_extent_rad
Definition: ERF_ReadBndryPlanes.H:96
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:90
@ RhoQ1_bc_comp
Definition: ERF_IndexDefines.H:91
@ RhoKE_bc_comp
Definition: ERF_IndexDefines.H:89
@ RhoTheta_bc_comp
Definition: ERF_IndexDefines.H:88
@ RhoQ2_bc_comp
Definition: ERF_IndexDefines.H:92
@ Rho_bc_comp
Definition: ERF_IndexDefines.H:87
@ xvel_bc
Definition: ERF_IndexDefines.H:102

Referenced by read_input_files().

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

◆ read_input_files()

void ReadBndryPlanes::read_input_files ( amrex::Real  time,
amrex::Real  dt,
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max m_bc_extdir_vals 
)

Function in ReadBndryPlanes for reading boundary data at a specific time and at the next timestep from input files.

Parameters
timeCurrent time
dtCurrent timestep
m_bc_extdir_valsContainer storing the external dirichlet boundary conditions we are reading from the input files
337 {
338  BL_PROFILE("ERF::ReadBndryPlanes::read_input_files");
339 
340  // Assert that both the current time and the next time are within the bounds
341  // of the data that we can read
342  AMREX_ALWAYS_ASSERT((m_in_times[0] <= time) && (time <= m_in_times.back()));
343  AMREX_ALWAYS_ASSERT((m_in_times[0] <= time+dt) && (time+dt <= m_in_times.back()));
344 
345  int ncomp = 1;
346 
347  const Box& domain = m_geom.Domain();
348  BoxArray ba(domain);
349  DistributionMapping dm{ba};
350  BndryRegister bndryn(ba, dm, m_in_rad, m_out_rad, m_extent_rad, ncomp);
351  bndryn.setVal(bogus_large_value);
352 
353  // The first time we enter this routine we read the first three files
354  if (last_file_read == -1)
355  {
356  int idx_init = 0;
357  read_file(idx_init,m_data_n ,m_bc_extdir_vals);
358  read_file(idx_init,m_data_interp,m_bc_extdir_vals); // We want to start with this filled
359  m_tn = m_in_times[idx_init];
360 
361  idx_init = 1;
362  read_file(idx_init,m_data_np1,m_bc_extdir_vals);
363  m_tnp1 = m_in_times[idx_init];
364 
365  idx_init = 2;
366  read_file(idx_init,m_data_np2,m_bc_extdir_vals);
367  m_tnp2 = m_in_times[idx_init];
368 
369  last_file_read = idx_init;
370  }
371 
372  // Compute the index such that time falls between times[idx] and times[idx+1]
373  const int idx = closest_index(m_in_times, time);
374 
375  // Advance the read window until it spans the requested time.
376  while (idx >= last_file_read-1 && last_file_read != m_in_times.size()-1) {
377  int new_read = last_file_read+1;
378 
379  // We need to change which data the pointers point to before we read in the new data
380  // This doesn't actually move the data, just swaps the pointers
381  for (OrientationIter oit; oit != nullptr; ++oit) {
382  auto ori = oit();
383  std::swap(m_data_n[ori] ,m_data_np1[ori]);
384  std::swap(m_data_np1[ori],m_data_np2[ori]);
385  }
386 
387  // Set the times corresponding to the post-swap pointers
388  m_tn = m_tnp1;
389  m_tnp1 = m_tnp2;
390  m_tnp2 = m_in_times[new_read];
391 
392  read_file(new_read,m_data_np2,m_bc_extdir_vals);
393  last_file_read = new_read;
394  }
395 
396  AMREX_ASSERT(time >= m_tn && time <= m_tnp2);
397  AMREX_ASSERT(time+dt >= m_tn && time+dt <= m_tnp2);
398 }
AMREX_FORCE_INLINE int closest_index(const Vector< Real > &vec, const Real value)
Definition: ERF_ReadBndryPlanes.cpp:15
amrex::Vector< amrex::Real > m_in_times
The timesteps / times that we read from time.dat.
Definition: ERF_ReadBndryPlanes.H:87
void read_file(int idx, amrex::Vector< std::unique_ptr< PlaneVector >> &data_to_fill, amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max > m_bc_extdir_vals)
Definition: ERF_ReadBndryPlanes.cpp:408
Here is the call graph for this function:

◆ read_time_file()

void ReadBndryPlanes::read_time_file ( )

Function in ReadBndryPlanes class for reading the external file specifying time data and broadcasting this data across MPI ranks.

265 {
266  BL_PROFILE("ERF::ReadBndryPlanes::read_time_file");
267 
268  // *********************************************************
269  // Read the time.data file and store the timesteps and times
270  // *********************************************************
271  int time_file_length = 0;
272 
273  if (ParallelDescriptor::IOProcessor()) {
274 
275  std::string line;
276  std::ifstream time_file(m_time_file);
277  if (!time_file.good()) {
278  Abort("Cannot find time file: " + m_time_file);
279  }
280  while (std::getline(time_file, line)) {
281  ++time_file_length;
282  }
283 
284  time_file.close();
285  }
286 
287  ParallelDescriptor::Bcast(
288  &time_file_length, 1,
289  ParallelDescriptor::IOProcessorNumber(),
290  ParallelDescriptor::Communicator());
291 
292  m_in_times.resize(time_file_length);
293  m_in_timesteps.resize(time_file_length);
294 
295  if (ParallelDescriptor::IOProcessor()) {
296  std::ifstream time_file(m_time_file);
297  for (int i = 0; i < time_file_length; ++i) {
298  time_file >> m_in_timesteps[i] >> m_in_times[i];
299  }
300  // Sanity check that there are no duplicates or mis-orderings
301  for (int i = 1; i < time_file_length; ++i) {
302  if (m_in_timesteps[i] <= m_in_timesteps[i-1])
303  Error("Bad timestep in time.dat file");
304  if (m_in_times[i] <= m_in_times[i-1])
305  Error("Bad time in time.dat file");
306  }
307  time_file.close();
308  }
309 
310  ParallelDescriptor::Bcast(
311  m_in_timesteps.data(), time_file_length,
312  ParallelDescriptor::IOProcessorNumber(),
313  ParallelDescriptor::Communicator());
314 
315  ParallelDescriptor::Bcast(
316  m_in_times.data(), time_file_length,
317  ParallelDescriptor::IOProcessorNumber(),
318  ParallelDescriptor::Communicator());
319 
320  // Allocate data we will need -- for now just at one level
321  int lev = 0;
322  define_level_data(lev);
323  Print() << "Successfully read time file and allocated data" << std::endl;
324 }
void define_level_data(int lev)
Definition: ERF_ReadBndryPlanes.cpp:43
Here is the call graph for this function:

◆ tinterp()

amrex::Real ReadBndryPlanes::tinterp ( ) const
inline
42 { return m_tinterp; }

Member Data Documentation

◆ is_density_read

int ReadBndryPlanes::is_density_read
private

◆ is_KE_read

int ReadBndryPlanes::is_KE_read
private

Referenced by ingested_KE(), and ReadBndryPlanes().

◆ is_q1_read

int ReadBndryPlanes::is_q1_read
private

Referenced by ingested_q1(), and ReadBndryPlanes().

◆ is_q2_read

int ReadBndryPlanes::is_q2_read
private

Referenced by ingested_q2(), and ReadBndryPlanes().

◆ is_scalar_read

int ReadBndryPlanes::is_scalar_read
private

Referenced by ingested_scalar(), and ReadBndryPlanes().

◆ is_temperature_read

int ReadBndryPlanes::is_temperature_read
private

Referenced by ingested_theta(), and ReadBndryPlanes().

◆ is_theta_read

int ReadBndryPlanes::is_theta_read
private

Referenced by ingested_theta(), and ReadBndryPlanes().

◆ is_velocity_read

int ReadBndryPlanes::is_velocity_read
private

◆ last_file_read

int ReadBndryPlanes::last_file_read
private

◆ m_data_interp

amrex::Vector<std::unique_ptr<PlaneVector> > ReadBndryPlanes::m_data_interp
private

Data interpolated to the time requested.

Referenced by read_input_files(), and ReadBndryPlanes().

◆ m_data_n

amrex::Vector<std::unique_ptr<PlaneVector> > ReadBndryPlanes::m_data_n
private

Data at time m_tn.

Referenced by read_input_files(), and ReadBndryPlanes().

◆ m_data_np1

amrex::Vector<std::unique_ptr<PlaneVector> > ReadBndryPlanes::m_data_np1
private

Data at time m_tnp1.

Referenced by read_input_files(), and ReadBndryPlanes().

◆ m_data_np2

amrex::Vector<std::unique_ptr<PlaneVector> > ReadBndryPlanes::m_data_np2
private

Data at time m_tnp2.

Referenced by read_input_files(), and ReadBndryPlanes().

◆ m_data_tendency

amrex::Vector<std::unique_ptr<PlaneVector> > ReadBndryPlanes::m_data_tendency
private

Tendency between the n and np1 data.

Referenced by ReadBndryPlanes().

◆ m_extent_rad

const int ReadBndryPlanes::m_extent_rad = 0
private

Referenced by read_file(), and read_input_files().

◆ m_filename

std::string ReadBndryPlanes::m_filename {""}
private

File name for IO.

Referenced by read_file(), and ReadBndryPlanes().

◆ m_geom

amrex::Geometry ReadBndryPlanes::m_geom
private

Geometry at level 0.

Referenced by read_file(), and read_input_files().

◆ m_in_rad

int ReadBndryPlanes::m_in_rad = 1
private

Controls extents on native bndry output.

Referenced by read_file(), read_input_files(), and ReadBndryPlanes().

◆ m_in_times

amrex::Vector<amrex::Real> ReadBndryPlanes::m_in_times
private

The timesteps / times that we read from time.dat.

Referenced by read_input_files(), and read_time_file().

◆ m_in_timesteps

amrex::Vector<int> ReadBndryPlanes::m_in_timesteps
private

Referenced by read_file(), and read_time_file().

◆ m_out_rad

const int ReadBndryPlanes::m_out_rad = 1
private

Referenced by read_file(), and read_input_files().

◆ m_rdOcp

const amrex::Real ReadBndryPlanes::m_rdOcp
private

R_d/c_p is needed for reading boundary files.

Referenced by read_file().

◆ m_time_file

std::string ReadBndryPlanes::m_time_file {""}
private

File name for file holding timesteps and times.

Referenced by read_time_file(), and ReadBndryPlanes().

◆ m_tinterp

amrex::Real ReadBndryPlanes::m_tinterp {-one}
private

Time for plane at interpolation.

Referenced by ReadBndryPlanes(), and tinterp().

◆ m_tn

amrex::Real ReadBndryPlanes::m_tn
private

The times for which we currently have data.

Referenced by read_input_files().

◆ m_tnp1

amrex::Real ReadBndryPlanes::m_tnp1
private

Referenced by read_input_files().

◆ m_tnp2

amrex::Real ReadBndryPlanes::m_tnp2
private

Referenced by read_input_files().

◆ m_use_real_bcs

bool ReadBndryPlanes::m_use_real_bcs = false
private

Are real BCs being used?

Referenced by read_file(), and ReadBndryPlanes().

◆ m_var_names

amrex::Vector<std::string> ReadBndryPlanes::m_var_names
private

Variables to be read in.

Referenced by read_file(), and ReadBndryPlanes().


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