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:7
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  const int t_step = m_in_timesteps[idx];
413  const std::string chkname1 = m_filename + Concatenate("/bndry_output", t_step);
414 
415  const std::string level_prefix = "Level_";
416  const int lev = 0;
417 
418  const Box& domain = m_geom.Domain();
419  BoxArray ba(domain);
420  DistributionMapping dm{ba};
421 
422  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NBCVAR_max> l_bc_extdir_vals_d;
423 
424  for (int i = 0; i < BCVars::NumTypes; i++)
425  {
426  for (OrientationIter oit; oit != nullptr; ++oit) {
427  auto ori = oit();
428  l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[i][ori];
429  }
430  }
431 
432  int n_for_density = -1;
433  for (int i = 0; i < m_var_names.size(); i++)
434  {
435  if (m_var_names[i] == "density") n_for_density = i;
436  }
437 
438  // We need to initialize all the components because we may not fill all of them from files,
439  // but the loop in the interpolate routine goes over all the components anyway
440  int ncomp_for_bc = BCVars::NumTypes;
441  for (OrientationIter oit; oit != nullptr; ++oit) {
442  auto ori = oit();
443  if (ori.coordDir() < 2) {
444  FArrayBox& d = (*data_to_fill[ori])[lev];
445  const auto& bx = d.box();
446  Array4<Real> d_arr = d.array();
447  ParallelFor(
448  bx, ncomp_for_bc, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
449  d_arr(i,j,k,n) = zero;
450  });
451  }
452  }
453 
454  // Read density for primitive to conserved conversions
455  std::string filenamer = MultiFabFileFullPrefix(lev, chkname1, level_prefix, "density");
456  BndryRegister bndry_r(ba, dm, m_in_rad, m_out_rad, m_extent_rad, 1);
457  bndry_r.setVal(Real(1.e13));
458  for (OrientationIter oit; oit != nullptr; ++oit) {
459  auto ori = oit();
460  if (ori.coordDir() < 2) {
461  std::string facenamer = Concatenate(filenamer + '_', ori, 1);
462  bndry_r[ori].read(facenamer);
463  }
464  }
465 
466  // Expose for GPU
467  bool real_bcs = m_use_real_bcs;
468 
469  for (int ivar = 0; ivar < m_var_names.size(); ivar++)
470  {
471  std::string var_name = m_var_names[ivar];
472 
473  std::string filename1 = MultiFabFileFullPrefix(lev, chkname1, level_prefix, var_name);
474 
475  int ncomp;
476  if (var_name == "velocity") {
477  ncomp = AMREX_SPACEDIM;
478  } else {
479  ncomp = 1;
480  }
481 
482  int n_offset;
483  if (var_name == "density") n_offset = BCVars::Rho_bc_comp;
484  if (var_name == "theta") n_offset = BCVars::RhoTheta_bc_comp;
485  if (var_name == "temperature") n_offset = BCVars::RhoTheta_bc_comp;
486  if (var_name == "ke") n_offset = BCVars::RhoKE_bc_comp;
487  if (var_name == "scalar") n_offset = BCVars::RhoScalar_bc_comp;
488  if (var_name == "qv") n_offset = BCVars::RhoQ1_bc_comp;
489  if (var_name == "qc") n_offset = BCVars::RhoQ2_bc_comp;
490  if (var_name == "velocity") n_offset = BCVars::xvel_bc;
491 
492  // Print() << "Reading " << chkname1 << " for variable " << var_name << " with n_offset == " << n_offset << std::endl;
493 
494  BndryRegister bndry(ba, dm, m_in_rad, m_out_rad, m_extent_rad, ncomp);
495  bndry.setVal(Real(1.e13));
496 
497  // *********************************************************
498  // Read in the BndryReg for all non-z faces
499  // *********************************************************
500  for (OrientationIter oit; oit != nullptr; ++oit) {
501  auto ori = oit();
502  if (ori.coordDir() < 2) {
503 
504  std::string facename1 = Concatenate(filename1 + '_', ori, 1);
505  bndry[ori].read(facename1);
506 
507  int normal = ori.coordDir();
508  IntVect v_offset = offset(ori.faceDir(), normal);
509  if (real_bcs) { v_offset = IntVect(0); }
510 
511  const auto& bbx = (*data_to_fill[ori])[lev].box();
512 
513  // *********************************************************
514  // Copy from the BndryReg into a MultiFab then use copyTo
515  // to write from the MultiFab to a single FAB for each face
516  // *********************************************************
517  MultiFab bndryMF(
518  bndry[ori].boxArray(), bndry[ori].DistributionMap(),
519  ncomp, 0, MFInfo());
520 
521  for (MFIter mfi(bndryMF); mfi.isValid(); ++mfi) {
522 
523  const auto& vbx = mfi.validbox();
524  const auto& bndry_read_arr = bndry[ori].array(mfi);
525  const auto& bndry_read_r_arr = bndry_r[ori].array(mfi);
526  const auto& bndry_mf_arr = bndryMF.array(mfi);
527 
528  const auto& bx = bbx & vbx;
529  if (bx.isEmpty()) {
530  continue;
531  }
532 
533  // We average the two cell-centered data points in the normal direction
534  // to define a Dirichlet value on the face itself.
535 
536  // This is the scalars -- they all get multiplied by rho, and in the case of
537  // reading in temperature, we must convert to theta first
538  Real rdOcp = m_rdOcp;
539  if (n_for_density >= 0) {
540  if (var_name == "temperature") {
541  ParallelFor(
542  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
543  Real R1 = bndry_read_r_arr(i, j, k, 0);
544  Real R2 = bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
545  Real T1 = bndry_read_arr(i, j, k, 0);
546  Real T2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
547  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
548  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
549  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
550  myhalf * (R1*Th1 + R2*Th2);
551  });
552  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
553  var_name == "qv" || var_name == "qc") {
554  ParallelFor(
555  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
556  Real R1 = bndry_read_r_arr(i, j, k, 0);
557  Real R2 = bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
558  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
559  myhalf * ( R1 * bndry_read_arr(i, j, k, 0) +
560  R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
561  });
562  } else if (var_name == "density") {
563  ParallelFor(
564  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
565  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
566  myhalf * ( bndry_read_arr(i, j, k, 0) +
567  bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
568  });
569  }
570  } else if (!ingested_density()) {
571  if (var_name == "temperature") {
572  ParallelFor(
573  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
574  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
575  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
576  Real T1 = bndry_read_arr(i, j, k, 0);
577  Real T2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0);
578  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
579  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
580  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
581  myhalf * (R1*Th1 + R2*Th2);
582  });
583  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
584  var_name == "qv" || var_name == "qc") {
585  ParallelFor(
586  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
587  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
588  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
589  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
590  myhalf * (R1 * bndry_read_arr(i, j, k, 0) +
591  R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
592  });
593  }
594  }
595 
596  // This is velocity
597  if (var_name == "velocity") {
598  ParallelFor(
599  bx, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
600  bndry_mf_arr(i, j, k, n) = (real_bcs) ? bndry_read_arr(i, j, k, n) :
601  myhalf * (bndry_read_arr(i, j, k, n) +
602  bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], n));
603  });
604  }
605 
606  } // mfi
607  bndryMF.copyTo((*data_to_fill[ori])[lev], 0, n_offset, ncomp);
608  } // coordDir < 2
609  } // ori
610  } // var_name
611 }
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
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(Real(1.e13));
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: