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::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::Real m_tinterp {-1.0}
 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$)
141 :
142  m_geom(geom),
143  m_rdOcp(rdOcp_in)
144 {
145  ParmParse pp("erf");
146 
147  // Get the radius inside the domain
148  pp.query("in_rad",m_in_rad);
149 
150  // Are we using real bcs?
151  pp.query("use_real_bcs", m_use_real_bcs);
152 
153  last_file_read = -1;
154 
155  m_tinterp = -1.;
156 
157  // What folder will the time series of planes be read from
158  pp.get("bndry_file", m_filename);
159 
160  is_velocity_read = 0;
161  is_density_read = 0;
163  is_theta_read = 0;
164  is_scalar_read = 0;
165  is_q1_read = 0;
166  is_q2_read = 0;
167  is_KE_read = 0;
168 
169  if (pp.contains("bndry_input_var_names"))
170  {
171  int num_vars = pp.countval("bndry_input_var_names");
172  m_var_names.resize(num_vars);
173  pp.queryarr("bndry_input_var_names",m_var_names,0,num_vars);
174  for (int i = 0; i < m_var_names.size(); i++) {
175  if (m_var_names[i] == "velocity") is_velocity_read = 1;
176  if (m_var_names[i] == "density") is_density_read = 1;
177  if (m_var_names[i] == "temperature") is_temperature_read = 1;
178  if (m_var_names[i] == "theta") is_theta_read = 1;
179  if (m_var_names[i] == "scalar") is_scalar_read = 1;
180  if (m_var_names[i] == "qv") is_q1_read = 1;
181  if (m_var_names[i] == "qc") is_q2_read = 1;
182  if (m_var_names[i] == "ke") is_KE_read = 1;
183  }
184  }
185 
186  // time.dat will be in the same folder as the time series of data
187  m_time_file = m_filename + "/time.dat";
188 
189  // each pointer (at at given time) has 6 components, one for each orientation
190  // TODO: we really only need 4 not 6
191  int size = 2*AMREX_SPACEDIM;
192  m_data_n.resize(size);
193  m_data_np1.resize(size);
194  m_data_np2.resize(size);
195  m_data_interp.resize(size);
196 }
ParmParse pp("prob")
int is_velocity_read
Definition: ERF_ReadBndryPlanes.H:98
int is_q2_read
Definition: ERF_ReadBndryPlanes.H:104
int is_theta_read
Definition: ERF_ReadBndryPlanes.H:101
bool m_use_real_bcs
Are real BCs being used?
Definition: ERF_ReadBndryPlanes.H:93
std::string m_filename
File name for IO.
Definition: ERF_ReadBndryPlanes.H:75
amrex::Real m_tinterp
Time for plane at interpolation.
Definition: ERF_ReadBndryPlanes.H:69
int is_temperature_read
Definition: ERF_ReadBndryPlanes.H:100
int is_density_read
Definition: ERF_ReadBndryPlanes.H:99
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_np2
Data at time m_tnp2.
Definition: ERF_ReadBndryPlanes.H:63
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_interp
Data interpolated to the time requested.
Definition: ERF_ReadBndryPlanes.H:66
int last_file_read
Definition: ERF_ReadBndryPlanes.H:107
int is_KE_read
Definition: ERF_ReadBndryPlanes.H:105
const amrex::Real m_rdOcp
R_d/c_p is needed for reading boundary files.
Definition: ERF_ReadBndryPlanes.H:96
std::string m_time_file
File name for file holding timesteps and times.
Definition: ERF_ReadBndryPlanes.H:78
amrex::Vector< std::string > m_var_names
Variables to be read in.
Definition: ERF_ReadBndryPlanes.H:85
int is_scalar_read
Definition: ERF_ReadBndryPlanes.H:102
int is_q1_read
Definition: ERF_ReadBndryPlanes.H:103
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_n
Data at time m_tn.
Definition: ERF_ReadBndryPlanes.H:57
amrex::Vector< std::unique_ptr< PlaneVector > > m_data_np1
Data at time m_tnp1.
Definition: ERF_ReadBndryPlanes.H:60
int m_in_rad
Controls extents on native bndry output.
Definition: ERF_ReadBndryPlanes.H:88
amrex::Geometry m_geom
Geometry at level 0.
Definition: ERF_ReadBndryPlanes.H:72
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 
60  const auto& lo = domain.loVect();
61  const auto& hi = domain.hiVect();
62 
63  IntVect plo(lo);
64  IntVect phi(hi);
65  const int normal = ori.coordDir();
66  plo[normal] = ori.isHigh() ? hi[normal] - (m_in_rad - 1) : -m_out_rad;
67  phi[normal] = ori.isHigh() ? hi[normal] + (m_out_rad ) : (m_in_rad - 1);
68  const Box pbx(plo, phi);
69  m_data_n[ori]->push_back(FArrayBox(pbx, ncomp));
70  m_data_np1[ori]->push_back(FArrayBox(pbx, ncomp));
71  m_data_np2[ori]->push_back(FArrayBox(pbx, ncomp));
72  m_data_interp[ori]->push_back(FArrayBox(pbx, ncomp));
73  }
74  }
75 }
const int m_out_rad
Definition: ERF_ReadBndryPlanes.H:89
@ NumTypes
Definition: ERF_IndexDefines.H:90

Referenced by read_time_file().

Here is the caller graph for this function:

◆ ingested_density()

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

Referenced by read_file().

Here is the caller graph for this function:

◆ ingested_KE()

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

◆ ingested_q1()

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

◆ ingested_q2()

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

◆ ingested_scalar()

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

◆ ingested_theta()

int ReadBndryPlanes::ingested_theta ( ) const
inline

◆ ingested_velocity()

int ReadBndryPlanes::ingested_velocity ( ) const
inline
41 {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
85 {
86  AMREX_ALWAYS_ASSERT(m_tn <= time && time <= m_tnp2);
87 
88  //Print() << "interp_in_time at time " << time << " given " << m_tn << " " << m_tnp1 << " " << m_tnp2 << std::endl;
89  //Print() << "m_tinterp " << m_tinterp << std::endl;
90 
91  if (time == m_tinterp) {
92  // We have already interpolated to this time
93  return m_data_interp;
94 
95  } else {
96 
97  // We must now interpolate to a new time
98  m_tinterp = time;
99 
100  if (time < m_tnp1) {
101  for (OrientationIter oit; oit != nullptr; ++oit) {
102  auto ori = oit();
103  if (ori.coordDir() < 2) {
104  const int nlevels = m_data_n[ori]->size();
105  for (int lev = 0; lev < nlevels; ++lev) {
106  const auto& datn = (*m_data_n[ori])[lev];
107  const auto& datnp1 = (*m_data_np1[ori])[lev];
108  auto& dati = (*m_data_interp[ori])[lev];
109  dati.linInterp<RunOn::Device>(
110  datn, 0, datnp1, 0, m_tn, m_tnp1, m_tinterp, datn.box(), 0, dati.nComp());
111  }
112  }
113  }
114  } else {
115  for (OrientationIter oit; oit != nullptr; ++oit) {
116  auto ori = oit();
117  if (ori.coordDir() < 2) {
118  const int nlevels = m_data_n[ori]->size();
119  for (int lev = 0; lev < nlevels; ++lev) {
120  const auto& datnp1 = (*m_data_np1[ori])[lev];
121  const auto& datnp2 = (*m_data_np2[ori])[lev];
122  auto& dati = (*m_data_interp[ori])[lev];
123  dati.linInterp<RunOn::Device>(
124  datnp1, 0, datnp2, 0, m_tnp1, m_tnp2, m_tinterp, datnp1.box(), 0,
125  dati.nComp());
126  }
127  }
128  }
129  }
130  }
131  return m_data_interp;
132 }
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
amrex::Real m_tnp1
Definition: ERF_ReadBndryPlanes.H:53
amrex::Real m_tnp2
Definition: ERF_ReadBndryPlanes.H:54
amrex::Real m_tn
The times for which we currently have data.
Definition: ERF_ReadBndryPlanes.H:52
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
349 {
350  const int t_step = m_in_timesteps[idx];
351  const std::string chkname1 = m_filename + Concatenate("/bndry_output", t_step);
352 
353  const std::string level_prefix = "Level_";
354  const int lev = 0;
355 
356  const Box& domain = m_geom.Domain();
357  BoxArray ba(domain);
358  DistributionMapping dm{ba};
359 
360  GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NBCVAR_max> l_bc_extdir_vals_d;
361 
362  for (int i = 0; i < BCVars::NumTypes; i++)
363  {
364  for (OrientationIter oit; oit != nullptr; ++oit) {
365  auto ori = oit();
366  l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[i][ori];
367  }
368  }
369 
370  int n_for_density = -1;
371  for (int i = 0; i < m_var_names.size(); i++)
372  {
373  if (m_var_names[i] == "density") n_for_density = i;
374  }
375 
376  // We need to initialize all the components because we may not fill all of them from files,
377  // but the loop in the interpolate routine goes over all the components anyway
378  int ncomp_for_bc = BCVars::NumTypes;
379  for (OrientationIter oit; oit != nullptr; ++oit) {
380  auto ori = oit();
381  if (ori.coordDir() < 2) {
382  FArrayBox& d = (*data_to_fill[ori])[lev];
383  const auto& bx = d.box();
384  Array4<Real> d_arr = d.array();
385  ParallelFor(
386  bx, ncomp_for_bc, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
387  d_arr(i,j,k,n) = 0.;
388  });
389  }
390  }
391 
392  // Read density for primitive to conserved conversions
393  std::string filenamer = MultiFabFileFullPrefix(lev, chkname1, level_prefix, "density");
394  BndryRegister bndry_r(ba, dm, m_in_rad, m_out_rad, m_extent_rad, 1);
395  bndry_r.setVal(1.0e13);
396  for (OrientationIter oit; oit != nullptr; ++oit) {
397  auto ori = oit();
398  if (ori.coordDir() < 2) {
399  std::string facenamer = Concatenate(filenamer + '_', ori, 1);
400  bndry_r[ori].read(facenamer);
401  }
402  }
403 
404  // Expose for GPU
405  bool real_bcs = m_use_real_bcs;
406 
407  for (int ivar = 0; ivar < m_var_names.size(); ivar++)
408  {
409  std::string var_name = m_var_names[ivar];
410 
411  std::string filename1 = MultiFabFileFullPrefix(lev, chkname1, level_prefix, var_name);
412 
413  int ncomp;
414  if (var_name == "velocity") {
415  ncomp = AMREX_SPACEDIM;
416  } else {
417  ncomp = 1;
418  }
419 
420  int n_offset;
421  if (var_name == "density") n_offset = BCVars::Rho_bc_comp;
422  if (var_name == "theta") n_offset = BCVars::RhoTheta_bc_comp;
423  if (var_name == "temperature") n_offset = BCVars::RhoTheta_bc_comp;
424  if (var_name == "ke") n_offset = BCVars::RhoKE_bc_comp;
425  if (var_name == "scalar") n_offset = BCVars::RhoScalar_bc_comp;
426  if (var_name == "qv") n_offset = BCVars::RhoQ1_bc_comp;
427  if (var_name == "qc") n_offset = BCVars::RhoQ2_bc_comp;
428  if (var_name == "velocity") n_offset = BCVars::xvel_bc;
429 
430  // Print() << "Reading " << chkname1 << " for variable " << var_name << " with n_offset == " << n_offset << std::endl;
431 
432  BndryRegister bndry(ba, dm, m_in_rad, m_out_rad, m_extent_rad, ncomp);
433  bndry.setVal(1.0e13);
434 
435  // *********************************************************
436  // Read in the BndryReg for all non-z faces
437  // *********************************************************
438  for (OrientationIter oit; oit != nullptr; ++oit) {
439  auto ori = oit();
440  if (ori.coordDir() < 2) {
441 
442  std::string facename1 = Concatenate(filename1 + '_', ori, 1);
443  bndry[ori].read(facename1);
444 
445  int normal = ori.coordDir();
446  IntVect v_offset = offset(ori.faceDir(), normal);
447  if (real_bcs) { v_offset = IntVect(0); }
448 
449  const auto& bbx = (*data_to_fill[ori])[lev].box();
450 
451  // *********************************************************
452  // Copy from the BndryReg into a MultiFab then use copyTo
453  // to write from the MultiFab to a single FAB for each face
454  // *********************************************************
455  MultiFab bndryMF(
456  bndry[ori].boxArray(), bndry[ori].DistributionMap(),
457  ncomp, 0, MFInfo());
458 
459  for (MFIter mfi(bndryMF); mfi.isValid(); ++mfi) {
460 
461  const auto& vbx = mfi.validbox();
462  const auto& bndry_read_arr = bndry[ori].array(mfi);
463  const auto& bndry_read_r_arr = bndry_r[ori].array(mfi);
464  const auto& bndry_mf_arr = bndryMF.array(mfi);
465 
466  const auto& bx = bbx & vbx;
467  if (bx.isEmpty()) {
468  continue;
469  }
470 
471  // We average the two cell-centered data points in the normal direction
472  // to define a Dirichlet value on the face itself.
473 
474  // This is the scalars -- they all get multiplied by rho, and in the case of
475  // reading in temperature, we must convert to theta first
476  Real rdOcp = m_rdOcp;
477  if (n_for_density >= 0) {
478  if (var_name == "temperature") {
479  ParallelFor(
480  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
481  Real R1 = bndry_read_r_arr(i, j, k, 0);
482  Real R2 = bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
483  Real T1 = bndry_read_arr(i, j, k, 0);
484  Real T2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
485  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
486  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
487  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
488  0.5 * (R1*Th1 + R2*Th2);
489  });
490  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
491  var_name == "qv" || var_name == "qc") {
492  ParallelFor(
493  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
494  Real R1 = bndry_read_r_arr(i, j, k, 0);
495  Real R2 = bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
496  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
497  0.5 * ( R1 * bndry_read_arr(i, j, k, 0) +
498  R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
499  });
500  } else if (var_name == "density") {
501  ParallelFor(
502  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
503  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
504  0.5 * ( bndry_read_arr(i, j, k, 0) +
505  bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
506  });
507  }
508  } else if (!ingested_density()) {
509  if (var_name == "temperature") {
510  ParallelFor(
511  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
512  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
513  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
514  Real T1 = bndry_read_arr(i, j, k, 0);
515  Real T2 = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0);
516  Real Th1 = getThgivenRandT(R1,T1,rdOcp);
517  Real Th2 = getThgivenRandT(R2,T2,rdOcp);
518  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
519  0.5 * (R1*Th1 + R2*Th2);
520  });
521  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
522  var_name == "qv" || var_name == "qc") {
523  ParallelFor(
524  bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
525  Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
526  Real R2 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
527  bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
528  0.5 * (R1 * bndry_read_arr(i, j, k, 0) +
529  R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
530  });
531  }
532  }
533 
534  // This is velocity
535  if (var_name == "velocity") {
536  ParallelFor(
537  bx, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
538  bndry_mf_arr(i, j, k, n) = (real_bcs) ? bndry_read_arr(i, j, k, n) :
539  0.5 * (bndry_read_arr(i, j, k, n) +
540  bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], n));
541  });
542  }
543 
544  } // mfi
545  bndryMF.copyTo((*data_to_fill[ori])[lev], 0, n_offset, ncomp);
546  } // coordDir < 2
547  } // ori
548  } // var_name
549 }
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=0.0)
Definition: ERF_EOS.H:64
#define NBCVAR_max
Definition: ERF_IndexDefines.H:29
const Real rdOcp
Definition: ERF_InitCustomPert_Bomex.H:16
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+0.5) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=1.0) { Real dT=T_pert *(std::cos(PI *L)+1.0)/2.0;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
amrex::Real Real
Definition: ERF_ShocInterface.H:19
amrex::Vector< int > m_in_timesteps
Definition: ERF_ReadBndryPlanes.H:82
int ingested_density() const
Definition: ERF_ReadBndryPlanes.H:43
const int m_extent_rad
Definition: ERF_ReadBndryPlanes.H:90
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:80
@ RhoQ1_bc_comp
Definition: ERF_IndexDefines.H:81
@ RhoKE_bc_comp
Definition: ERF_IndexDefines.H:79
@ RhoTheta_bc_comp
Definition: ERF_IndexDefines.H:78
@ RhoQ2_bc_comp
Definition: ERF_IndexDefines.H:82
@ Rho_bc_comp
Definition: ERF_IndexDefines.H:77
@ xvel_bc
Definition: ERF_IndexDefines.H:87

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
275 {
276  BL_PROFILE("ERF::ReadBndryPlanes::read_input_files");
277 
278  // Assert that both the current time and the next time are within the bounds
279  // of the data that we can read
280  AMREX_ALWAYS_ASSERT((m_in_times[0] <= time) && (time <= m_in_times.back()));
281  AMREX_ALWAYS_ASSERT((m_in_times[0] <= time+dt) && (time+dt <= m_in_times.back()));
282 
283  int ncomp = 1;
284 
285  const Box& domain = m_geom.Domain();
286  BoxArray ba(domain);
287  DistributionMapping dm{ba};
288  BndryRegister bndryn(ba, dm, m_in_rad, m_out_rad, m_extent_rad, ncomp);
289  bndryn.setVal(1.0e13);
290 
291  // The first time we enter this routine we read the first three files
292  if (last_file_read == -1)
293  {
294  int idx_init = 0;
295  read_file(idx_init,m_data_n,m_bc_extdir_vals);
296  read_file(idx_init,m_data_interp,m_bc_extdir_vals); // We want to start with this filled
297  m_tn = m_in_times[idx_init];
298 
299  idx_init = 1;
300  read_file(idx_init,m_data_np1,m_bc_extdir_vals);
301  m_tnp1 = m_in_times[idx_init];
302 
303  idx_init = 2;
304  read_file(idx_init,m_data_np2,m_bc_extdir_vals);
305  m_tnp2 = m_in_times[idx_init];
306 
307  last_file_read = idx_init;
308  }
309 
310  // Compute the index such that time falls between times[idx] and times[idx+1]
311  const int idx = closest_index(m_in_times, time);
312 
313  // Advance the read window until it spans the requested time.
314  while (idx >= last_file_read-1 && last_file_read != m_in_times.size()-1) {
315  int new_read = last_file_read+1;
316 
317  // We need to change which data the pointers point to before we read in the new data
318  // This doesn't actually move the data, just swaps the pointers
319  for (OrientationIter oit; oit != nullptr; ++oit) {
320  auto ori = oit();
321  std::swap(m_data_n[ori] ,m_data_np1[ori]);
322  std::swap(m_data_np1[ori],m_data_np2[ori]);
323  }
324 
325  // Set the times corresponding to the post-swap pointers
326  m_tn = m_tnp1;
327  m_tnp1 = m_tnp2;
328  m_tnp2 = m_in_times[new_read];
329 
330  read_file(new_read,m_data_np2,m_bc_extdir_vals);
331  last_file_read = new_read;
332  }
333 
334  AMREX_ASSERT(time >= m_tn && time <= m_tnp2);
335  AMREX_ASSERT(time+dt >= m_tn && time+dt <= m_tnp2);
336 }
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:81
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:346
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.

203 {
204  BL_PROFILE("ERF::ReadBndryPlanes::read_time_file");
205 
206  // *********************************************************
207  // Read the time.data file and store the timesteps and times
208  // *********************************************************
209  int time_file_length = 0;
210 
211  if (ParallelDescriptor::IOProcessor()) {
212 
213  std::string line;
214  std::ifstream time_file(m_time_file);
215  if (!time_file.good()) {
216  Abort("Cannot find time file: " + m_time_file);
217  }
218  while (std::getline(time_file, line)) {
219  ++time_file_length;
220  }
221 
222  time_file.close();
223  }
224 
225  ParallelDescriptor::Bcast(
226  &time_file_length, 1,
227  ParallelDescriptor::IOProcessorNumber(),
228  ParallelDescriptor::Communicator());
229 
230  m_in_times.resize(time_file_length);
231  m_in_timesteps.resize(time_file_length);
232 
233  if (ParallelDescriptor::IOProcessor()) {
234  std::ifstream time_file(m_time_file);
235  for (int i = 0; i < time_file_length; ++i) {
236  time_file >> m_in_timesteps[i] >> m_in_times[i];
237  }
238  // Sanity check that there are no duplicates or mis-orderings
239  for (int i = 1; i < time_file_length; ++i) {
240  if (m_in_timesteps[i] <= m_in_timesteps[i-1])
241  Error("Bad timestep in time.dat file");
242  if (m_in_times[i] <= m_in_times[i-1])
243  Error("Bad time in time.dat file");
244  }
245  time_file.close();
246  }
247 
248  ParallelDescriptor::Bcast(
249  m_in_timesteps.data(), time_file_length,
250  ParallelDescriptor::IOProcessorNumber(),
251  ParallelDescriptor::Communicator());
252 
253  ParallelDescriptor::Bcast(
254  m_in_times.data(), time_file_length,
255  ParallelDescriptor::IOProcessorNumber(),
256  ParallelDescriptor::Communicator());
257 
258  // Allocate data we will need -- for now just at one level
259  int lev = 0;
260  define_level_data(lev);
261  Print() << "Successfully read time file and allocated data" << std::endl;
262 }
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
39 { 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_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 {-1.0}
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: