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 (double time, double 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 double &time)
amrex::Vector< std::unique_ptr< PlaneVector > > &	get_tendency (const double &time)
double	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

Static Public Member Functions
static int	plane_slot (const std::string &var_name)

Private Attributes
double	m_tn
	The times for which we currently have data. More...
double	m_tnp1
double	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...
double	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< double >	m_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

geom	Geometry for the domain
rdOcp_in	Real constant for the Rhydberg constant ($R_d$) divided by the specific heat at constant pressure ($c_p$)

:
    m_geom(geom),
    m_rdOcp(rdOcp_in)
{
    ParmParse pp("erf");
 
    // Get the radius inside the domain
    pp.queryAdd("in_rad",m_in_rad);
 
    // Are we using real bcs?
    pp.queryAdd("use_real_bcs", m_use_real_bcs);
 
    last_file_read = -1;
 
    m_tinterp = -1.0;
 
    // What folder will the time series of planes be read from
    pp.get("bndry_file", m_filename);
 
    is_velocity_read     = 0;
    is_density_read      = 0;
    is_temperature_read  = 0;
    is_theta_read        = 0;
    is_scalar_read       = 0;
    is_q1_read           = 0;
    is_q2_read           = 0;
    is_KE_read           = 0;
 
    if (pp.contains("bndry_input_var_names"))
    {
        int num_vars = pp.countval("bndry_input_var_names");
        m_var_names.resize(num_vars);
        pp.queryarr("bndry_input_var_names",m_var_names,0,num_vars);
        for (int i = 0; i < m_var_names.size(); i++) {
            if (m_var_names[i] == "velocity")     is_velocity_read = 1;
            if (m_var_names[i] == "density")      is_density_read = 1;
            if (m_var_names[i] == "temperature")  is_temperature_read = 1;
            if (m_var_names[i] == "theta")        is_theta_read = 1;
            if (m_var_names[i] == "scalar")       is_scalar_read = 1;
            if (m_var_names[i] == "qv")           is_q1_read = 1;
            if (m_var_names[i] == "qc")           is_q2_read = 1;
            if (m_var_names[i] == "ke")           is_KE_read = 1;
            if (plane_slot(m_var_names[i]) < 0) {
                Abort("erf.bndry_input_var_names: unknown variable \"" + m_var_names[i] +
                      "\"; the names read are velocity, density, temperature, theta, scalar, qv, qc and ke");
            }
        }
        // The passive scalars share one boundary condition, so an ingested scalar
        // plane fills every one of them; the writer saves only the first.
        // (if constexpr: with one scalar the body is compiled out, so clang's
        // -Wunreachable-code does not flag it.)
        if constexpr (NSCALARS > 1) {
            if (is_scalar_read) {
                Print() << "WARNING: erf.bndry_input_var_names reads \"scalar\": the plane of the first "
                        << "passive scalar (rhoadv_0) fills all " << NSCALARS << " passive scalars at the "
                        << "ingested faces, the extra ones (such as dust or smoke) included" << std::endl;
            }
        }
    }
 
    // time.dat will be in the same folder as the time series of data
    m_time_file = m_filename + "/time.dat";
 
    // each pointer (at at given time) has 6 components, one for each orientation
    // TODO: we really only need 4 not 6
    int size = 2*AMREX_SPACEDIM;
    m_data_n.resize(size);
    m_data_np1.resize(size);
    m_data_np2.resize(size);
    m_data_interp.resize(size);
    m_data_tendency.resize(size);
}

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.

{
    Print() << "ReadBndryPlanes::define_level_data" << std::endl;
    // *********************************************************
    // Allocate space for all of the boundary planes we may need
    // *********************************************************
    int ncomp = BCVars::NumTypes;
    const Box& domain = m_geom.Domain();
    for (OrientationIter oit; oit != nullptr; ++oit) {
        auto ori = oit();
        if (ori.coordDir() < 2) {
 
            m_data_n[ori]        = std::make_unique<PlaneVector>();
            m_data_np1[ori]      = std::make_unique<PlaneVector>();
            m_data_np2[ori]      = std::make_unique<PlaneVector>();
            m_data_interp[ori]   = std::make_unique<PlaneVector>();
            m_data_tendency[ori] = std::make_unique<PlaneVector>();
 
            const auto& lo = domain.loVect();
            const auto& hi = domain.hiVect();
 
            IntVect plo(lo);
            IntVect phi(hi);
            const int normal = ori.coordDir();
            plo[normal] = ori.isHigh() ? hi[normal] - (m_in_rad - 1) : lo[normal] - m_out_rad;
            phi[normal] = ori.isHigh() ? hi[normal] + (m_out_rad   ) : lo[normal] + (m_in_rad - 1);
            const Box pbx(plo, phi);
            m_data_n[ori]->push_back(FArrayBox(pbx, ncomp));
            m_data_np1[ori]->push_back(FArrayBox(pbx, ncomp));
            m_data_np2[ori]->push_back(FArrayBox(pbx, ncomp));
            m_data_interp[ori]->push_back(FArrayBox(pbx, ncomp));
            m_data_tendency[ori]->push_back(FArrayBox(pbx, ncomp));
        }
    }
}

Referenced by read_time_file().

Here is the call graph for this function:

Here is the caller graph for this function:

get_tendency()

Vector< std::unique_ptr< PlaneVector > > & ReadBndryPlanes::get_tendency

(

const double &

time_in

)

Function in ReadBndryPlanes class for interpolating boundary data in time.

Parameters

time_in

Constant specifying the time for interpolation

{
    // A restart that lands exactly on a boundary-plane time can request a time a few
    // ULP outside [m_tn, m_tnp2] because per-level t_new drifts under subcycling.
    // Tolerate that drift, then clamp into the valid window before interpolating.
    const double eps  = 1.0e-8 * (m_tnp2 - m_tn);
    AMREX_ALWAYS_ASSERT(m_tn - eps <= time_in && time_in <= m_tnp2 + eps);
    const double time = std::min(std::max(time_in, m_tn), m_tnp2);
 
    if (time < m_tnp1) {
        Real idt = static_cast<Real>(1.0 / (m_tnp1 - m_tn));
        for (OrientationIter oit; oit != nullptr; ++oit) {
            auto ori = oit();
            if (ori.coordDir() < 2) {
                const int nlevels = static_cast<int>(m_data_n[ori]->size());
                for (int lev = 0; lev < nlevels; ++lev) {
                    auto& fabt = (*m_data_tendency[ori])[lev];
                    Box bx     = fabt.box();
                    int ncomp  = fabt.nComp();
 
                    const auto& datt   = fabt.array();
                    const auto& datn   = (*m_data_n[ori])[lev].array();
                    const auto& datnp1 = (*m_data_np1[ori])[lev].array();
                    ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                    {
                        datt(i,j,k,n) = (datnp1(i,j,k,n) - datn(i,j,k,n)) * idt;
                    });
                }
            }
        }
    } else {
        Real idt = static_cast<Real>(1.0 / (m_tnp2 - m_tnp1));
        for (OrientationIter oit; oit != nullptr; ++oit) {
            auto ori = oit();
            if (ori.coordDir() < 2) {
                const int nlevels = static_cast<int>(m_data_n[ori]->size());
                for (int lev = 0; lev < nlevels; ++lev) {
                    auto& fabt = (*m_data_tendency[ori])[lev];
                    Box bx     = fabt.box();
                    int ncomp  = fabt.nComp();
 
                    const auto& datt   = fabt.array();
                    const auto& datnp1 = (*m_data_np1[ori])[lev].array();
                    const auto& datnp2 = (*m_data_np2[ori])[lev].array();
                    ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                    {
                        datt(i,j,k,n) = (datnp2(i,j,k,n) - datnp1(i,j,k,n)) * idt;
                    });
                }
            }
        }
    }
 
    return m_data_tendency;
}

Here is the call graph for this function:

ingested_density()

int ReadBndryPlanes::ingested_density ( ) const

inline

{return is_density_read;}

Referenced by read_file().

Here is the caller graph for this function:

ingested_KE()

int ReadBndryPlanes::ingested_KE ( ) const

inline

{return is_KE_read;}

ingested_q1()

int ReadBndryPlanes::ingested_q1 ( ) const

inline

{return is_q1_read;}

ingested_q2()

int ReadBndryPlanes::ingested_q2 ( ) const

inline

{return is_q2_read;}

ingested_scalar()

int ReadBndryPlanes::ingested_scalar ( ) const

inline

{return is_scalar_read;}

ingested_theta()

int ReadBndryPlanes::ingested_theta ( ) const

inline

{return (is_temperature_read || is_theta_read);}

ingested_velocity()

int ReadBndryPlanes::ingested_velocity ( ) const

inline

{return is_velocity_read;}

interp_in_time()

Vector< std::unique_ptr< PlaneVector > > & ReadBndryPlanes::interp_in_time

(

const double &

time_in

)

Function in ReadBndryPlanes class for interpolating boundary data in time.

Parameters

time_in

Constant specifying the time for interpolation

{
    // A restart that lands exactly on a boundary-plane time can request a time a few
    // ULP outside [m_tn, m_tnp2] because per-level t_new drifts under subcycling.
    // Tolerate that drift, then clamp into the valid window before interpolating.
    const double eps  = 1.0e-8 * (m_tnp2 - m_tn);
    AMREX_ALWAYS_ASSERT(m_tn - eps <= time_in && time_in <= m_tnp2 + eps);
    const double time = std::min(std::max(time_in, m_tn), m_tnp2);
 
    //Print() << "interp_in_time at time " << time << " given " << m_tn << " " << m_tnp1 << " " << m_tnp2 << std::endl;
    //Print() << "m_tinterp " << m_tinterp << std::endl;
 
    if (time == m_tinterp) {
        // We have already interpolated to this time
        return m_data_interp;
 
    } else {
 
        // We must now interpolate to a new time
        m_tinterp = time;
 
        if (time < m_tnp1) {
            for (OrientationIter oit; oit != nullptr; ++oit) {
                auto ori = oit();
                if (ori.coordDir() < 2) {
                    const int nlevels = static_cast<int>(m_data_n[ori]->size());
                    for (int lev = 0; lev < nlevels; ++lev) {
                        const auto& datn   = (*m_data_n[ori])[lev];
                        const auto& datnp1 = (*m_data_np1[ori])[lev];
                        auto& dati = (*m_data_interp[ori])[lev];
                        dati.linInterp<RunOn::Device>(datn, 0, datnp1, 0,
                                                      static_cast<Real>(m_tn), static_cast<Real>(m_tnp1),
                                                      static_cast<Real>(m_tinterp),
                                                      datn.box(), 0, dati.nComp());
                    }
                }
            }
        } else {
            for (OrientationIter oit; oit != nullptr; ++oit) {
                auto ori = oit();
                if (ori.coordDir() < 2) {
                    const int nlevels = static_cast<int>(m_data_n[ori]->size());
                    for (int lev = 0; lev < nlevels; ++lev) {
                        const auto& datnp1 = (*m_data_np1[ori])[lev];
                        const auto& datnp2 = (*m_data_np2[ori])[lev];
                        auto& dati = (*m_data_interp[ori])[lev];
                        dati.linInterp<RunOn::Device>(datnp1, 0, datnp2, 0,
                                                      static_cast<Real>(m_tnp1), static_cast<Real>(m_tnp2),
                                                      static_cast<Real>(m_tinterp),
                                                      datnp1.box(), 0, dati.nComp());
                    }
                }
            }
        }
    }
    return m_data_interp;
}

Here is the call graph for this function:

plane_slot()

static int ReadBndryPlanes::plane_slot ( const std::string &  var_name )

inlinestatic

The slot (BCVars) a plane of variable var_name is stored in, the same slot fill_from_bndryregs reads it from (cons_bc_comp of the component it fills); -1 for a name that is not a plane variable.

    {
        if (var_name == "density")     { return BCVars::Rho_bc_comp; }
        if (var_name == "theta")       { return BCVars::RhoTheta_bc_comp; }
        if (var_name == "temperature") { return BCVars::RhoTheta_bc_comp; }
        if (var_name == "ke")          { return BCVars::RhoKE_bc_comp; }
        if (var_name == "scalar")      { return BCVars::RhoScalar_bc_comp; }
        if (var_name == "qv")          { return BCVars::RhoQ1_bc_comp; }
        if (var_name == "qc")          { return BCVars::RhoQ2_bc_comp; }
        if (var_name == "velocity")    { return BCVars::xvel_bc; }
        return -1;
    }

Referenced by read_file(), and ReadBndryPlanes().

Here is the caller 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

idx	Specifies the index corresponding to the timestep we want
data_to_fill	Container for face data on boundaries
m_bc_extdir_vals	Container storing the external dirichlet boundary conditions we are reading from the input files

{
    if (idx >= m_in_timesteps.size()) {
        Print() << "Asking for index " << idx << " but m_in_timesteps only has size " << m_in_timesteps.size() << std::endl;
        Abort();
    }
    const int t_step = m_in_timesteps[idx];
    const std::string chkname1 = m_filename + Concatenate("/bndry_output", t_step);
 
    const std::string level_prefix = "Level_";
    const int lev = 0;
 
    const Box& domain = m_geom.Domain();
    BoxArray ba(domain);
    DistributionMapping dm{ba};
 
    GpuArray<GpuArray<Real, AMREX_SPACEDIM*2>, AMREX_SPACEDIM+NBCVAR_max> l_bc_extdir_vals_d;
 
    for (int i = 0; i < BCVars::NumTypes; i++)
    {
        for (OrientationIter oit; oit != nullptr; ++oit) {
            auto ori = oit();
            l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[i][ori];
        }
    }
 
    int n_for_density = -1;
    for (int i = 0; i < m_var_names.size(); i++)
    {
       if (m_var_names[i] == "density") n_for_density = i;
    }
 
    // We need to initialize all the components because we may not fill all of them from files,
    //    but the loop in the interpolate routine goes over all the components anyway
    int ncomp_for_bc = BCVars::NumTypes;
    for (OrientationIter oit; oit != nullptr; ++oit) {
        auto ori = oit();
        if (ori.coordDir() < 2) {
            FArrayBox& d = (*data_to_fill[ori])[lev];
            const auto& bx = d.box();
            Array4<Real> d_arr = d.array();
            ParallelFor(
                bx, ncomp_for_bc, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
                d_arr(i,j,k,n) = zero;
            });
        }
    }
 
    // Read density for primitive to conserved conversions
    std::string filenamer = MultiFabFileFullPrefix(lev, chkname1, level_prefix, "density");
    BndryRegister bndry_r(ba, dm, m_in_rad, m_out_rad, m_extent_rad, 1);
    bndry_r.setVal(bogus_large_value);
    for (OrientationIter oit; oit != nullptr; ++oit) {
          auto ori = oit();
          if (ori.coordDir() < 2) {
              std::string facenamer = Concatenate(filenamer + '_', ori, 1);
              bndry_r[ori].read(facenamer);
          }
    }
 
    // A moist writer computes the temperature plane with qv (erf_dermoisttemp), so the conversion
    // back to theta below needs the same qv: read its plane too when both are read.  Otherwise qv
    // is zero there, as for a dry writer.
    BndryRegister bndry_qv(ba, dm, m_in_rad, m_out_rad, m_extent_rad, 1);
    bndry_qv.setVal(zero);
    if (is_temperature_read && is_q1_read && !m_use_real_bcs) {
        std::string filenameqv = MultiFabFileFullPrefix(lev, chkname1, level_prefix, "qv");
        for (OrientationIter oit; oit != nullptr; ++oit) {
              auto ori = oit();
              if (ori.coordDir() < 2) {
                  std::string facenameqv = Concatenate(filenameqv + '_', ori, 1);
                  bndry_qv[ori].read(facenameqv);
              }
        }
    }
 
    // A temperature plane written with qv cannot be converted back to theta without it
    if (is_temperature_read && !is_q1_read && !m_use_real_bcs) {
        const std::string qv_header = Concatenate(MultiFabFileFullPrefix(lev, chkname1, level_prefix, "qv") + '_',
                                                  Orientation(0, Orientation::low), 1) + "_H";
        if (FileExists(qv_header)) {
            Abort("erf.bndry_input_var_names reads temperature but not qv, and " + m_filename +
                  " holds qv planes: the writer was moist, so its temperature includes qv.  Add qv to"
                  " erf.bndry_input_var_names, or write and read theta instead of temperature");
        }
    }
 
    // Expose for GPU
    bool real_bcs = m_use_real_bcs;
 
    for (int ivar = 0; ivar < m_var_names.size(); ivar++)
    {
        std::string var_name = m_var_names[ivar];
 
        std::string filename1 = MultiFabFileFullPrefix(lev, chkname1, level_prefix, var_name);
 
        int ncomp;
        if (var_name == "velocity") {
            ncomp = AMREX_SPACEDIM;
        } else {
            ncomp = 1;
        }
 
        // The slot the plane is stored in; every name was checked against
        // plane_slot in the constructor, so it is never -1 here.
        const int n_offset = plane_slot(var_name);
 
        // Print() << "Reading " << chkname1 << " for variable " << var_name << " with n_offset == " << n_offset << std::endl;
 
        BndryRegister bndry(ba, dm, m_in_rad, m_out_rad, m_extent_rad, ncomp);
        bndry.setVal(bogus_large_value);
 
        // *********************************************************
        // Read in the BndryReg for all non-z faces
        // *********************************************************
        for (OrientationIter oit; oit != nullptr; ++oit) {
          auto ori = oit();
          if (ori.coordDir() < 2) {
 
            std::string facename1 = Concatenate(filename1 + '_', ori, 1);
            bndry[ori].read(facename1);
 
            int normal = ori.coordDir();
            IntVect v_offset = offset(ori.faceDir(), normal);
            if (real_bcs) { v_offset = IntVect(0); }
 
            const auto& bbx = (*data_to_fill[ori])[lev].box();
 
            // *********************************************************
            // Copy from the BndryReg into a MultiFab then use copyTo
            //     to write from the MultiFab to a single FAB for each face
            // *********************************************************
            MultiFab bndryMF(
                bndry[ori].boxArray(), bndry[ori].DistributionMap(),
                ncomp, 0, MFInfo());
 
            for (MFIter mfi(bndryMF); mfi.isValid(); ++mfi) {
 
                const auto& vbx = mfi.validbox();
                const auto& bndry_read_arr   = bndry[ori].array(mfi);
                const auto& bndry_read_r_arr = bndry_r[ori].array(mfi);
                const auto& bndry_read_qv_arr = bndry_qv[ori].array(mfi);
                const auto& bndry_mf_arr     = bndryMF.array(mfi);
 
                const auto& bx = bbx & vbx;
                if (bx.isEmpty()) {
                    continue;
                }
 
                // Split the 2-cell-thick working box into ghost and interior
                // slots so the (i+v_offset) neighbor access stays in-bounds.
                // Both slots are filled with the same face-averaged Dirichlet
                // value; the interior slot just needs the opposite neighbor.
                Box bx_ghost = bx;
                Box bx_int   = bx;
                if (ori.isLow()) {
                    bx_ghost.setBig  (normal, domain.smallEnd(normal) - 1);
                    bx_int  .setSmall(normal, domain.smallEnd(normal));
                } else {
                    bx_ghost.setSmall(normal, domain.bigEnd(normal) + 1);
                    bx_int  .setBig  (normal, domain.bigEnd(normal));
                }
                const IntVect v_offset_int = -v_offset;
 
                // We average the two cell-centered data points in the normal direction
                //    to define a Dirichlet value on the face itself.
 
                // This is the scalars -- they all get multiplied by rho, and in the case of
                //   reading in temperature, we must convert to theta first
                Real rdOcp = m_rdOcp;
                if (n_for_density >= 0) {
                  if (var_name == "temperature") {
                    ParallelFor(
                        bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                             Real R1 =  bndry_read_r_arr(i, j, k, 0);
                             Real R2 =  bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
                             Real T1 =  bndry_read_arr(i, j, k, 0);
                             Real T2 =  bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
                             Real Q1  = bndry_read_qv_arr(i, j, k, 0);
                             Real Q2  = bndry_read_qv_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
                             Real Th1 = getThgivenRandT(R1,T1,rdOcp,Q1);
                             Real Th2 = getThgivenRandT(R2,T2,rdOcp,Q2);
                             bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                                                     myhalf * (R1*Th1 + R2*Th2);
                        });
                  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
                             var_name == "qv"    || var_name == "qc") {
                    ParallelFor(
                        bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                             Real R1 =  bndry_read_r_arr(i, j, k, 0);
                             Real R2 =  bndry_read_r_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
                             bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                 myhalf * ( R1 * bndry_read_arr(i, j, k, 0) +
                                         R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
                        });
                   } else if (var_name == "density") {
                    ParallelFor(
                        bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                                bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                    myhalf * ( bndry_read_arr(i, j, k, 0) +
                                            bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
                        });
                   }
                } else if (!ingested_density()) {
                  if (var_name == "temperature") {
                    ParallelFor(
                        bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                             Real R1  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                             Real R2  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                             Real T1  = bndry_read_arr(i, j, k, 0);
                             Real T2  = bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0);
                             Real Q1  = bndry_read_qv_arr(i, j, k, 0);
                             Real Q2  = bndry_read_qv_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2],0);
                             Real Th1 = getThgivenRandT(R1,T1,rdOcp,Q1);
                             Real Th2 = getThgivenRandT(R2,T2,rdOcp,Q2);
                             bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                                                     myhalf * (R1*Th1 + R2*Th2);
                        });
                  } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
                             var_name == "qv"    || var_name == "qc") {
                      ParallelFor(
                        bx_ghost, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                             Real R1  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                             Real R2  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                             bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                myhalf * (R1 * bndry_read_arr(i, j, k, 0) +
                                       R2 * bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], 0));
                        });
                  }
                }
 
                // This is velocity
                if (var_name == "velocity") {
                    ParallelFor(
                        bx_ghost, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
                                bndry_mf_arr(i, j, k, n) = (real_bcs) ? bndry_read_arr(i, j, k, n) :
                                  myhalf * (bndry_read_arr(i, j, k, n) +
                                         bndry_read_arr(i+v_offset[0],j+v_offset[1],k+v_offset[2], n));
                        });
                }
 
                // --- interior-slot fill (same face-averaged Dirichlet value;
                //     neighbor offset is flipped because the "other cell"
                //     is now on the opposite side of the boundary face) ---
                if (n_for_density >= 0) {
                    if (var_name == "temperature") {
                        ParallelFor(
                            bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                                 Real R1 =  bndry_read_r_arr(i, j, k, 0);
                                 Real R2 =  bndry_read_r_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
                                 Real T1 =  bndry_read_arr(i, j, k, 0);
                                 Real T2 =  bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
                                 Real Q1  = bndry_read_qv_arr(i, j, k, 0);
                                 Real Q2  = bndry_read_qv_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
                                 Real Th1 = getThgivenRandT(R1,T1,rdOcp,Q1);
                                 Real Th2 = getThgivenRandT(R2,T2,rdOcp,Q2);
                                 bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                                                         myhalf * (R1*Th1 + R2*Th2);
                            });
                    } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
                               var_name == "qv"    || var_name == "qc") {
                        ParallelFor(
                            bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                                 Real R1 =  bndry_read_r_arr(i, j, k, 0);
                                 Real R2 =  bndry_read_r_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
                                 bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                     myhalf * ( R1 * bndry_read_arr(i, j, k, 0) +
                                             R2 * bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0));
                            });
                    } else if (var_name == "density") {
                        ParallelFor(
                            bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                                    bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                        myhalf * ( bndry_read_arr(i, j, k, 0) +
                                                bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0));
                            });
                    }
                } else if (!ingested_density()) {
                    if (var_name == "temperature") {
                        ParallelFor(
                            bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                                 Real R1  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                                 Real R2  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                                 Real T1  = bndry_read_arr(i, j, k, 0);
                                 Real T2  = bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0);
                                 Real Q1  = bndry_read_qv_arr(i, j, k, 0);
                                 Real Q2  = bndry_read_qv_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2],0);
                                 Real Th1 = getThgivenRandT(R1,T1,rdOcp,Q1);
                                 Real Th2 = getThgivenRandT(R2,T2,rdOcp,Q2);
                                 bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                                                         myhalf * (R1*Th1 + R2*Th2);
                            });
                    } else if (var_name == "theta" || var_name == "ke" || var_name == "scalar" ||
                               var_name == "qv"    || var_name == "qc") {
                          ParallelFor(
                            bx_int, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                                 Real R1  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                                 Real R2  = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori];
                                 bndry_mf_arr(i, j, k, 0) = (real_bcs) ? bndry_read_arr(i, j, k, 0) :
                                    myhalf * (R1 * bndry_read_arr(i, j, k, 0) +
                                           R2 * bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], 0));
                            });
                    }
                }
 
                if (var_name == "velocity") {
                    ParallelFor(
                        bx_int, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
                                bndry_mf_arr(i, j, k, n) = (real_bcs) ? bndry_read_arr(i, j, k, n) :
                                  myhalf * (bndry_read_arr(i, j, k, n) +
                                         bndry_read_arr(i+v_offset_int[0],j+v_offset_int[1],k+v_offset_int[2], n));
                        });
                }
 
            } // mfi
            bndryMF.copyTo((*data_to_fill[ori])[lev], 0, n_offset, ncomp);
          } // coordDir < 2
        } // ori
    } // var_name
}

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	(	double	time,
		double	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

time	Current time
dt	Current timestep
m_bc_extdir_vals	Container storing the external dirichlet boundary conditions we are reading from the input files

{
    BL_PROFILE("ERF::ReadBndryPlanes::read_input_files");
 
    // Assert that both the current time and the next time are within the bounds
    // of the data that we can read
    AMREX_ALWAYS_ASSERT((m_in_times[0] <= time) && (time <= m_in_times.back()));
    AMREX_ALWAYS_ASSERT((m_in_times[0] <= time+dt) && (time+dt <= m_in_times.back()));
 
    int ncomp = 1;
 
    const Box& domain = m_geom.Domain();
    BoxArray ba(domain);
    DistributionMapping dm{ba};
    BndryRegister bndryn(ba, dm, m_in_rad, m_out_rad, m_extent_rad, ncomp);
    bndryn.setVal(bogus_large_value);
 
    // The first time we enter this routine we read the first three files
    if (last_file_read == -1)
    {
        int idx_init = 0;
        read_file(idx_init,m_data_n     ,m_bc_extdir_vals);
        read_file(idx_init,m_data_interp,m_bc_extdir_vals); // We want to start with this filled
        m_tn = m_in_times[idx_init];
 
        idx_init = 1;
        read_file(idx_init,m_data_np1,m_bc_extdir_vals);
        m_tnp1 = m_in_times[idx_init];
 
        idx_init = 2;
        read_file(idx_init,m_data_np2,m_bc_extdir_vals);
        m_tnp2 = m_in_times[idx_init];
 
        last_file_read = idx_init;
    }
 
    // Compute the index such that time falls between times[idx] and times[idx+1]
    const int idx = closest_index(m_in_times, time);
 
    // Advance the read window until it spans the requested time.
    while (idx >= last_file_read-1 && last_file_read != m_in_times.size()-1) {
        int new_read = last_file_read+1;
 
        // We need to change which data the pointers point to before we read in the new data
        // This doesn't actually move the data, just swaps the pointers
        for (OrientationIter oit; oit != nullptr; ++oit) {
            auto ori = oit();
            std::swap(m_data_n[ori]  ,m_data_np1[ori]);
            std::swap(m_data_np1[ori],m_data_np2[ori]);
        }
 
        // Set the times corresponding to the post-swap pointers
        m_tn   = m_tnp1;
        m_tnp1 = m_tnp2;
        m_tnp2 = m_in_times[new_read];
 
        read_file(new_read,m_data_np2,m_bc_extdir_vals);
        last_file_read = new_read;
    }
 
    AMREX_ASSERT(time    >= m_tn && time    <= m_tnp2);
    AMREX_ASSERT(time+dt >= m_tn && time+dt <= m_tnp2);
}

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.

{
    BL_PROFILE("ERF::ReadBndryPlanes::read_time_file");
 
    // *********************************************************
    // Read the time.data file and store the timesteps and times
    // *********************************************************
    int time_file_length = 0;
 
    if (ParallelDescriptor::IOProcessor()) {
 
        std::string line;
        std::ifstream time_file(m_time_file);
        if (!time_file.good()) {
            Abort("Cannot find time file: " + m_time_file);
        }
        while (std::getline(time_file, line)) {
            ++time_file_length;
        }
 
        time_file.close();
    }
 
    ParallelDescriptor::Bcast(
        &time_file_length, 1,
        ParallelDescriptor::IOProcessorNumber(),
        ParallelDescriptor::Communicator());
 
    m_in_times.resize(time_file_length);
    m_in_timesteps.resize(time_file_length);
 
    if (ParallelDescriptor::IOProcessor()) {
        std::ifstream time_file(m_time_file);
        for (int i = 0; i < time_file_length; ++i) {
            time_file >> m_in_timesteps[i] >> m_in_times[i];
        }
        // Sanity check that there are no duplicates or mis-orderings
        for (int i = 1; i < time_file_length; ++i) {
            if (m_in_timesteps[i] <= m_in_timesteps[i-1])
                Error("Bad timestep in time.dat file");
            if (m_in_times[i] <= m_in_times[i-1])
                Error("Bad time in time.dat file");
        }
        time_file.close();
    }
 
    ParallelDescriptor::Bcast(
        m_in_timesteps.data(), time_file_length,
        ParallelDescriptor::IOProcessorNumber(),
        ParallelDescriptor::Communicator());
 
    ParallelDescriptor::Bcast(
        m_in_times.data(), time_file_length,
        ParallelDescriptor::IOProcessorNumber(),
        ParallelDescriptor::Communicator());
 
    // Allocate data we will need -- for now just at one level
    int lev = 0;
    define_level_data(lev);
    Print() << "Successfully read time file and allocated data" << std::endl;
}

Here is the call graph for this function:

tinterp()

double ReadBndryPlanes::tinterp ( ) const

inline

{ return m_tinterp; }

Member Data Documentation

is_density_read

int ReadBndryPlanes::is_density_read

private

Referenced by ingested_density(), and ReadBndryPlanes().

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(), read_file(), 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(), read_file(), 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

Referenced by ingested_velocity(), and ReadBndryPlanes().

last_file_read

int ReadBndryPlanes::last_file_read

private

Referenced by read_input_files(), and ReadBndryPlanes().

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<double> 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

double ReadBndryPlanes::m_tinterp {-1.0}

private

Time for plane at interpolation.

Referenced by ReadBndryPlanes(), and tinterp().

m_tn

double ReadBndryPlanes::m_tn

private

The times for which we currently have data.

Referenced by read_input_files().

m_tnp1

double ReadBndryPlanes::m_tnp1

private

Referenced by read_input_files().

m_tnp2

double 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:

• Source/IO/ERF_ReadBndryPlanes.H
• Source/IO/ERF_ReadBndryPlanes.cpp