#include <ERF_SampleData.H>

Collaboration diagram for PlaneSampler:

[legend]

Public Member Functions
	PlaneSampler ()

amrex::Box	getIndexBox (const amrex::RealBox &real_box, const amrex::Geometry &geom)

void	get_sample_data (amrex::Vector< amrex::Geometry > &geom, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_new)

void	write_sample_data (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::RealBox >	m_bnd_rbx

amrex::Vector< std::unique_ptr< amrex::MultiFab > >	m_ps_mf

amrex::Vector< std::string >	m_name

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

Constructor & Destructor Documentation

◆ PlaneSampler()

PlaneSampler::PlaneSampler ( )

inline

     {
         amrex::ParmParse pp("erf");
  
         // Count number of lo and hi points define the plane
         int n_plane_lo  = pp.countval("sample_plane_lo") / AMREX_SPACEDIM;
         int n_plane_hi  = pp.countval("sample_plane_hi") / AMREX_SPACEDIM;
         int n_plane_dir = pp.countval("sample_plane_dir");
         AMREX_ALWAYS_ASSERT( (n_plane_lo==n_plane_hi ) &&
                              (n_plane_lo==n_plane_dir) );
  
         // Parse the data
         if (n_plane_lo > 0) {
             // Parse lo
             amrex::Vector<amrex::Real> r_lo; r_lo.resize(n_plane_lo*AMREX_SPACEDIM);
             amrex::Vector<amrex::Vector<amrex::Real>> rv_lo;
             pp.queryarr("sample_plane_lo",r_lo,0,n_plane_lo*AMREX_SPACEDIM);
             for (int i(0); i < n_plane_lo; i++) {
                 amrex::Vector<amrex::Real> rv = {r_lo[AMREX_SPACEDIM*i+0],
                                                  r_lo[AMREX_SPACEDIM*i+1],
                                                  r_lo[AMREX_SPACEDIM*i+2]};
                 rv_lo.push_back(rv);
             }
  
             // Parse hi
             amrex::Vector<amrex::Real> r_hi; r_hi.resize(n_plane_hi*AMREX_SPACEDIM);
             amrex::Vector<amrex::Vector<amrex::Real>> rv_hi;
             pp.queryarr("sample_plane_hi",r_hi,0,n_plane_hi*AMREX_SPACEDIM);
             for (int i(0); i < n_plane_hi; i++) {
                 amrex::Vector<amrex::Real> rv = {r_hi[AMREX_SPACEDIM*i+0],
                                                  r_hi[AMREX_SPACEDIM*i+1],
                                                  r_hi[AMREX_SPACEDIM*i+2]};
                 rv_hi.push_back(rv);
             }
  
             // Construct vector of bounding real boxes
             m_bnd_rbx.resize(n_plane_lo);
             for (int i(0); i < n_plane_hi; i++){
                 amrex::RealBox rbx(rv_lo[i].data(),rv_hi[i].data());
                 m_bnd_rbx[i] = rbx;
             }
  
             // Parse directionality
             m_dir.resize(n_plane_dir);
             pp.queryarr("sample_plane_dir",m_dir,0,n_plane_dir);
  
             // Parse names
             std::string name_base = "plt_plane_";
             m_name.resize(n_plane_lo);
             int n_names = pp.countval("sample_plane_name");
             if (n_names > 0) {
                 AMREX_ALWAYS_ASSERT( n_names==n_plane_lo );
                 pp.queryarr("sample_plane_name",m_name,0,n_names);
             } else {
                 for (int iplane(0); iplane<n_plane_lo; ++iplane) {
                     m_name[iplane] = amrex::Concatenate(name_base, iplane , 5);
                 }
             }
  
             // Allocate space for level indicator
             m_lev.resize(n_plane_dir,0);
  
             // Allocate space for MF pointers
             m_ps_mf.resize(n_plane_lo);
  
             // Get requested vars
             if (pp.countval("plane_sampling_vars") > 0) {
                 m_varnames.clear();
                 amrex::Vector<std::string> requested_vars;
                 pp.queryarr("plane_sampling_vars",requested_vars);
                 amrex::Print() << "Selected plane sampling vars :";
                 if (containerHasElement(requested_vars, "density")) {
                     m_varnames.push_back("density");
                     amrex::Print() << " " << "density";
                 }
                 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, "qc")) {
                     m_varnames.push_back("qc");
                     amrex::Print() << " " << "qc";
                 }
                 if (containerHasElement(requested_vars, "pressure")) {
                     m_varnames.push_back("pressure");
                     amrex::Print() << " " << "pressure";
                 }
                 amrex::Print() << std::endl;
             }
         }
     }

Here is the call graph for this function:

Member Function Documentation

◆ get_sample_data()

void PlaneSampler::get_sample_data	(	amrex::Vector< amrex::Geometry > &	geom,
		amrex::Vector< amrex::Vector< amrex::MultiFab >> &	vars_new
	)

inline

     {
         int nlev   = static_cast<int>(vars_new.size());
         int nplane = static_cast<int>(m_bnd_rbx.size());
         int ncomp  = static_cast<int>(m_varnames.size());
         bool interpolate = true;
  
         int qv_comp = -1;
  
         // Loop over each plane
         for (int iplane(0); iplane<nplane; ++iplane) {
             int dir = m_dir[iplane];
             amrex::RealBox bnd_rbx  = m_bnd_rbx[iplane];
             amrex::Real point = bnd_rbx.lo(dir);
  
             // 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 MultiFab holding requested variables
                 amrex::MultiFab mf_cc_data;
                 mf_cc_data.define(ba, dm, ncomp, 1);
  
                 int mf_comp = 0;
  
                 if (containerHasElement(m_varnames, "density")) {
                     amrex::MultiFab::Copy(mf_cc_data, vars_new[ilev][Vars::cons], Rho_comp, mf_comp, 1, 0);
                     mf_comp += 1;
                 }
  
                 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")) {
                     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp,
                         "qv sampling requested but moisture components not present in state");
                     if (qv_comp >= 0) AMREX_ALWAYS_ASSERT(qv_comp == mf_comp);
                     qv_comp = mf_comp;
 #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;
                 }
                 if (containerHasElement(m_varnames, "qc")) {
                     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ2_comp,
                         "qc sampling requested but moisture components not present in state");
 #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,RhoQ2_comp) / cfab(i,j,k,Rho_comp);
                         });
                     }
                     mf_comp += 1;
                 }
  
                 if (containerHasElement(m_varnames, "pressure")) {
                     if (vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp) {
                         // moist pressure: use qv from dfab if already sampled, else compute inline
 #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
                             {
                                 amrex::Real qv_val = (qv_comp >= 0) ? dfab(i,j,k,qv_comp)
                                                                      : cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
                                 dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp), qv_val);
                             });
                         }
                     } else {
                         // dry pressure
 #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[iplane] = ilev;
                 m_ps_mf[iplane] = get_slice_data(dir, point, mf_cc_data, geom[ilev],
                                                  0, ncomp, interpolate, bnd_rbx);
  
                 // We can stop if we got the entire plane
                 auto min_bnd_bx   = m_ps_mf[iplane]->boxArray().minimalBox();
                 amrex::Box bnd_bx = getIndexBox(bnd_rbx, geom[ilev]);
                 if (bnd_bx == min_bnd_bx) { break; }
  
             } // ilev
         }// iplane
     }

Here is the call graph for this function:

◆ getIndexBox()

amrex::Box PlaneSampler::getIndexBox	(	const amrex::RealBox &	real_box,
		const amrex::Geometry &	geom
	)

inline

                                             {
         amrex::IntVect slice_lo, slice_hi;
  
         AMREX_D_TERM(slice_lo[0]=static_cast<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
                      slice_lo[1]=static_cast<int>(std::floor((real_box.lo(1) - geom.ProbLo(1))/geom.CellSize(1)));,
                      slice_lo[2]=static_cast<int>(std::floor((real_box.lo(2) - geom.ProbLo(2))/geom.CellSize(2))););
  
         AMREX_D_TERM(slice_hi[0]=static_cast<int>(std::floor((real_box.hi(0) - geom.ProbLo(0))/geom.CellSize(0)));,
                      slice_hi[1]=static_cast<int>(std::floor((real_box.hi(1) - geom.ProbLo(1))/geom.CellSize(1)));,
                      slice_hi[2]=static_cast<int>(std::floor((real_box.hi(2) - geom.ProbLo(2))/geom.CellSize(2))););
  
         return amrex::Box(slice_lo, slice_hi) & geom.Domain();
     }

Referenced by get_sample_data(), and write_sample_data().

Here is the caller graph for this function:

◆ write_sample_data()

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

inline

     {
         int nplane = m_ps_mf.size();
         for (int iplane(0); iplane<nplane; ++iplane) {
             // Data members that can be used as-is
             int dir = m_dir[iplane];
             int lev = m_lev[iplane];
             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 plane
             amrex::RealBox m_rb = m_bnd_rbx[iplane];
             amrex::Box m_dom    = getIndexBox(m_rb, geom[lev]);
             amrex::Real point   = m_rb.hi(dir);
             amrex::Vector<int> is_per(AMREX_SPACEDIM,0);
             for (int d(0); d<AMREX_SPACEDIM; ++d) {
                 if (d==dir) {
                     m_rb.setLo(d,point-myhalf*geom[lev].CellSize(d));
                     m_rb.setHi(d,point+myhalf*geom[lev].CellSize(d));
                 }
                 is_per[d] = geom[lev].isPeriodic(d);
             }
             amrex::Vector<amrex::Geometry> m_geom; m_geom.resize(1);
             m_geom[0].define(m_dom, &m_rb, geom[lev].Coord(), is_per.data());
  
             // Create plotfile name
             std::string name_plane = m_name[iplane];
             name_plane += "_step_";
             std::string plotfilename = amrex::Concatenate(name_plane, m_level_steps[0], 5);
  
             // Get the data
             amrex::Vector<const amrex::MultiFab*> mf = {m_ps_mf[iplane].get()};
  
             // Write each plane
             WriteMultiLevelPlotfile(plotfilename, 1, mf,
                                     m_varnames, m_geom, m_time,
                                     m_level_steps, m_ref_ratio);
         } // iplane
     }