#include <ERF_SampleData.H>

Collaboration diagram for LineSampler:

[legend]

Public Member Functions
	LineSampler ()

void	get_line_mfs (amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_new)

void	write_line_mfs (amrex::Vector< amrex::Real > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)

void	write_line_ascii (amrex::Vector< amrex::Real > &time)

void	write_line_plotfile (amrex::Vector< amrex::Real > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)

Public Attributes
amrex::Vector< int >	m_dir

amrex::Vector< int >	m_lev

amrex::Vector< amrex::Box >	m_bnd_bx

amrex::Vector< amrex::MultiFab >	m_ls_mf

amrex::Vector< std::string >	m_name

bool	m_write_ascii {false}

amrex::Vector< std::string >	m_varnames {"magvel","theta"}

amrex::Vector< std::unique_ptr< std::fstream > >	m_datastream

Constructor & Destructor Documentation

◆ LineSampler()

LineSampler::LineSampler ( )

inline

     {
         amrex::ParmParse pp("erf");
  
         // Count number of lo and hi points define the line
         int n_line_lo  = pp.countval("sample_line_lo") / AMREX_SPACEDIM;
         int n_line_hi  = pp.countval("sample_line_hi") / AMREX_SPACEDIM;
         int n_line_dir = pp.countval("sample_line_dir");
         AMREX_ALWAYS_ASSERT( (n_line_lo==n_line_hi ) &&
                              (n_line_lo==n_line_dir) );
  
         // Parse the data
         if (n_line_lo > 0) {
             // Parse lo
             amrex::Vector<int> idx_lo; idx_lo.resize(n_line_lo*AMREX_SPACEDIM);
             amrex::Vector<amrex::IntVect> iv_lo; iv_lo.resize(n_line_lo);
             pp.queryarr("sample_line_lo",idx_lo,0,n_line_lo*AMREX_SPACEDIM);
             for (int i(0); i < n_line_lo; i++) {
                 amrex::IntVect iv(idx_lo[AMREX_SPACEDIM*i+0],
                                   idx_lo[AMREX_SPACEDIM*i+1],
                                   idx_lo[AMREX_SPACEDIM*i+2]);
                 iv_lo[i] = iv;
             }
  
             // Parse hi
             amrex::Vector<int> idx_hi; idx_hi.resize(n_line_hi*AMREX_SPACEDIM);
             amrex::Vector<amrex::IntVect> iv_hi; iv_hi.resize(n_line_hi);
             pp.queryarr("sample_line_hi",idx_hi,0,n_line_hi*AMREX_SPACEDIM);
             for (int i(0); i < n_line_hi; i++) {
                 amrex::IntVect iv(idx_hi[AMREX_SPACEDIM*i+0],
                                   idx_hi[AMREX_SPACEDIM*i+1],
                                   idx_hi[AMREX_SPACEDIM*i+2]);
                 iv_hi[i] = iv;
             }
  
             // Construct vector of bounding boxes
             m_bnd_bx.resize(n_line_lo);
             for (int i = 0; i < n_line_hi; i++){
                 amrex::Box lbx(iv_lo[i],iv_hi[i]);
                 m_bnd_bx[i] = lbx;
             }
  
             // Parse directionality
             m_dir.resize(n_line_dir);
             pp.queryarr("sample_line_dir",m_dir,0,n_line_dir);
  
             // Parse names
             std::string name_base = "plt_line_";
             m_name.resize(n_line_lo);
             int n_names = pp.countval("sample_line_name");
             if (n_names > 0) {
                 AMREX_ALWAYS_ASSERT( n_names==n_line_lo );
                 pp.queryarr("sample_line_name",m_name,0,n_names);
             } else {
                 for (int iline(0); iline<n_line_lo; ++iline) {
                     m_name[iline] = amrex::Concatenate(name_base, iline , 5);
                 }
             }
  
             // Allocate space for level indicator
             m_lev.resize(n_line_dir,0);
  
             // Allocate space for MF pointers
             m_ls_mf.resize(n_line_lo);
  
             // Get requested vars
             if (pp.countval("line_sampling_vars") > 0) {
                 m_varnames.clear();
                 amrex::Vector<std::string> requested_vars;
                 pp.queryarr("line_sampling_vars",requested_vars);
                 amrex::Print() << "Selected line sampling vars :";
                 if (containerHasElement(requested_vars, "x_velocity")) {
                     m_varnames.push_back("x_velocity");
                     amrex::Print() << " " << "x_velocity";
                 }
                 if (containerHasElement(requested_vars, "y_velocity")) {
                     m_varnames.push_back("y_velocity");
                     amrex::Print() << " " << "y_velocity";
                 }
                 if (containerHasElement(requested_vars, "z_velocity")) {
                     m_varnames.push_back("z_velocity");
                     amrex::Print() << " " << "z_velocity";
                 }
                 if (containerHasElement(requested_vars, "magvel")) {
                     m_varnames.push_back("magvel");
                     amrex::Print() << " " << "magvel";
                 }
                 if (containerHasElement(requested_vars, "theta")) {
                     m_varnames.push_back("theta");
                     amrex::Print() << " " << "theta";
                 }
                 if (containerHasElement(requested_vars, "qv")) {
                     m_varnames.push_back("qv");
                     amrex::Print() << " " << "qv";
                 }
                 if (containerHasElement(requested_vars, "pressure")) {
                     m_varnames.push_back("pressure");
                     amrex::Print() << " " << "pressure";
                 }
                 amrex::Print() << std::endl;
             }
  
             // Write outputs to text files, one file per variable, with all
             // times appended to the same file
             pp.query("line_sampling_text_output",m_write_ascii);
             if (m_write_ascii && amrex::ParallelDescriptor::IOProcessor()) {
                 int nvar = static_cast<int>(m_varnames.size());
                 m_datastream.resize(n_line_lo * nvar);
                 int i = 0;
                 for (int iline(0); iline<n_line_lo; ++iline) {
                     for (int ivar(0); ivar<nvar; ++ivar) {
                         std::string filename = m_name[iline] + "." + m_varnames[ivar];
                         m_datastream[i] = std::make_unique<std::fstream>();
                         m_datastream[i]->open(filename.c_str(),std::ios::out|std::ios::app);
                         if (!m_datastream[i]->good()) {
                             amrex::FileOpenFailed(filename);
                         }
                         i++;
                     }
                 }
             }
         }
     }

Here is the call graph for this function:

Member Function Documentation

◆ get_line_mfs()

void LineSampler::get_line_mfs ( amrex::Vector< amrex::Vector< amrex::MultiFab >> & vars_new )

inline

     {
         int nlev  = static_cast<int>(vars_new.size());
         int nline = static_cast<int>(m_bnd_bx.size());
         int ncomp = static_cast<int>(m_varnames.size());
  
         // Loop over each line
         for (int iline(0); iline<nline; ++iline) {
             int dir = m_dir[iline];
             amrex::Box bnd_bx   = m_bnd_bx[iline];
             amrex::IntVect cell = bnd_bx.smallEnd();
  
             // Search each level to get the finest data possible
             for (int ilev(nlev-1); ilev>=0; --ilev) {
  
                 // Construct CC velocities
                 amrex::MultiFab mf_cc_vel;
                 auto ba = vars_new[ilev][Vars::cons].boxArray();
                 auto dm = vars_new[ilev][Vars::cons].DistributionMap();
                 mf_cc_vel.define(ba, dm, AMREX_SPACEDIM, amrex::IntVect(1,1,1));
                 average_face_to_cellcenter(mf_cc_vel,0,
                                            amrex::Array<const amrex::MultiFab*,3>{&vars_new[ilev][Vars::xvel],
                                                                                   &vars_new[ilev][Vars::yvel],
                                                                                   &vars_new[ilev][Vars::zvel]});
  
                 // Construct vector of MFs holding T and WSP
                 amrex::MultiFab mf_cc_data;
                 mf_cc_data.define(ba, dm, ncomp, 1);
  
                 int mf_comp = 0;
  
                 if (containerHasElement(m_varnames, "x_velocity")) {
                     amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 0, mf_comp, 1, 0);
                     mf_comp += 1;
                 }
                 if (containerHasElement(m_varnames, "y_velocity")) {
                     amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 1, mf_comp, 1, 0);
                     mf_comp += 1;
                 }
                 if (containerHasElement(m_varnames, "z_velocity")) {
                     amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 2, mf_comp, 1, 0);
                     mf_comp += 1;
                 }
  
                 if (containerHasElement(m_varnames, "magvel")) {
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
                     for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
                         const amrex::Box& tbx = mfi.tilebox();
                         auto const& dfab = mf_cc_data.array(mfi);
                         auto const& vfab = mf_cc_vel.array(mfi);
  
                         amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                         {
                             dfab(i,j,k,mf_comp) = std::sqrt(vfab(i,j,k,0)*vfab(i,j,k,0)
                                                           + vfab(i,j,k,1)*vfab(i,j,k,1)
                                                           + vfab(i,j,k,2)*vfab(i,j,k,2)) ;
                         });
                     }
                     mf_comp += 1;
                 }
  
                 if (containerHasElement(m_varnames, "theta")) {
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
                     for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
                         const amrex::Box& tbx = mfi.tilebox();
                         auto const& dfab = mf_cc_data.array(mfi);
                         auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
  
                         amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                         {
                             dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoTheta_comp) / cfab(i,j,k,Rho_comp);
                         });
                     }
                     mf_comp += 1;
                 }
  
                 if (containerHasElement(m_varnames, "qv") && containerHasElement(m_varnames, "pressure")) {
                     // if qv is requested, assume that we have moisture
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
                     for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
                         const amrex::Box& tbx = mfi.tilebox();
                         auto const& dfab = mf_cc_data.array(mfi);
                         auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
  
                         amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                         {
                             dfab(i,j,k,mf_comp  ) = cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
                             dfab(i,j,k,mf_comp+1) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp),
                                                                  dfab(i,j,k,mf_comp));
                         });
                     }
                     mf_comp += 2;
                 } else if (containerHasElement(m_varnames, "qv")) {
                     // if qv is requested, assume that we have moisture
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
                     for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
                         const amrex::Box& tbx = mfi.tilebox();
                         auto const& dfab = mf_cc_data.array(mfi);
                         auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
  
                         amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                         {
                             dfab(i,j,k,mf_comp  ) = cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
                         });
                     }
                     mf_comp += 1;
                 } else if (containerHasElement(m_varnames, "pressure")) {
                     // pressure requested w/o qv, assume dry
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
                     for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
                         const amrex::Box& tbx = mfi.tilebox();
                         auto const& dfab = mf_cc_data.array(mfi);
                         auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
  
                         amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
                         {
                             dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp));
                         });
                     }
                     mf_comp += 1;
                 }
  
                   m_lev[iline] = ilev;
                 m_ls_mf[iline] = get_line_data(mf_cc_data, dir, cell, bnd_bx);
  
                 // We can stop if we got the entire line
                 auto min_bnd_bx = m_ls_mf[iline].boxArray().minimalBox();
                 if (bnd_bx == min_bnd_bx) { continue; }
  
             } // ilev
         }// iline
     }

Here is the call graph for this function:

◆ write_line_ascii()

void LineSampler::write_line_ascii ( amrex::Vector< amrex::Real > & time )

inline

     {
         // same as definitions in ERF.H
         constexpr int datwidth = 14;
         constexpr int datprecision = 6;
         constexpr int timeprecision = 13;
  
         int nline = static_cast<int>(m_ls_mf.size());
         int nvar  = static_cast<int>(m_varnames.size());
         for (int iline(0); iline<nline; ++iline) {
             int dir = m_dir[iline];
             int lev = m_lev[iline];
             amrex::Real m_time = time[lev];
             amrex::Box m_dom = m_bnd_bx[iline];
  
             for (int ivar(0); ivar<nvar; ++ivar) {
                 // Convert multifab to vector
                 amrex::Gpu::HostVector<amrex::Real> vec = sumToLine(m_ls_mf[iline], ivar, 1, m_dom, dir);
  
                 if (amrex::ParallelDescriptor::IOProcessor()) {
                     int ifile = iline*nvar + ivar;
                     std::ostream& fs = *m_datastream[ifile];
                     fs << std::setw(datwidth) << std::setprecision(timeprecision) << m_time
                        << std::setw(datwidth) << std::setprecision(datprecision);
                     for (const auto& val : vec) {
                         fs << " " << val;
                     }
                     fs << std::endl;
                 }
             }
         }
     }

Referenced by write_line_mfs().

Here is the caller graph for this function:

◆ write_line_mfs()

void LineSampler::write_line_mfs	(	amrex::Vector< amrex::Real > &	time,
		amrex::Vector< int > &	level_steps,
		amrex::Vector< amrex::IntVect > &	ref_ratio,
		amrex::Vector< amrex::Geometry > &	geom
	)

inline

     {
         if (m_write_ascii) {
             write_line_ascii(time);
         } else {
             write_line_plotfile(time, level_steps, ref_ratio, geom);
         }
     }

Here is the call graph for this function:

◆ write_line_plotfile()

void LineSampler::write_line_plotfile	(	amrex::Vector< amrex::Real > &	time,
		amrex::Vector< int > &	level_steps,
		amrex::Vector< amrex::IntVect > &	ref_ratio,
		amrex::Vector< amrex::Geometry > &	geom
	)

inline

     {
         int nline = static_cast<int>(m_ls_mf.size());
         for (int iline(0); iline<nline; ++iline) {
             // Data members that can be used as-is
             int dir = m_dir[iline];
             int lev = m_lev[iline];
             amrex::Real m_time = time[lev];
             amrex::Vector<int> m_level_steps = {level_steps[lev]};
             amrex::Vector<amrex::IntVect> m_ref_ratio  = {ref_ratio[lev]};
  
             // Create modified geometry object corresponding to the line
             auto plo = geom[lev].ProbLo();
             auto dx  = geom[lev].CellSize();
             amrex::Vector<amrex::Geometry> m_geom; m_geom.resize(1);
             amrex::Vector<int> is_per(AMREX_SPACEDIM,0);
             amrex::Box m_dom = m_bnd_bx[iline];
             amrex::RealBox m_rb;
             for (int d(0); d<AMREX_SPACEDIM; ++d) {
                 amrex::Real offset = (d==dir) ? 0 : 0.5;
                 amrex::Real lo = plo[d] + ( m_dom.smallEnd(d) - offset ) * dx[d];
                 amrex::Real hi = plo[d] + ( m_dom.bigEnd(d)   + offset ) * dx[d];
  
                 m_rb.setLo(d,lo);
                 m_rb.setHi(d,hi);
  
                 is_per[d] = geom[lev].isPeriodic(d);
             }
             m_geom[0].define(m_dom, &m_rb, geom[lev].Coord(), is_per.data());
  
             // Create plotfile name
             std::string name_line = m_name[iline];
             name_line += "_step_";
             std::string plotfilename = amrex::Concatenate(name_line, m_level_steps[0], 5);
  
             // Get the data
             amrex::Vector<const amrex::MultiFab*> mf = {&(m_ls_mf[iline])};
  
             // Write each line
             WriteMultiLevelPlotfile(plotfilename, 1, mf,
                                     m_varnames, m_geom, m_time,
                                     m_level_steps, m_ref_ratio);
         }
     }