Super rough but it plotted a line.

Exec/ABL/inputs_sample

@@ -20,15 +20,19 @@ erf.fixed_dt           = 0.1  # fixed time step depending on grid resolution
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
 erf.v              = 1       # verbosity in ERF.cpp
-amr.v                = 1       # verbosity in Amr.cpp
+amr.v              = 1       # verbosity in Amr.cpp
 
 erf.data_log = my_data_file
 
 erf.sample_point_log = my_sample_point
 erf.sample_point = 6 63 5
 
 erf.sample_line_log = my_sample_line
-erf.sample_line = 63 63 5
+erf.sample_line = 5 32 5
+
+erf.sample_line_lo  = 5 32  5 10 32  5
+erf.sample_line_hi  = 5 32 25 10 32 25
+erf.sample_line_dir = 2 2
 
 # REFINEMENT / REGRIDDING
 amr.max_level       = 0       # maximum level number allowed

Source/ERF.H

@@ -40,6 +40,7 @@
 #include <ERF_MRI.H>
 #include <ERF_PhysBCFunct.H>
 #include <ERF_FillPatcher.H>
+#include <ERF_SampleData.H>
 
 #ifdef ERF_USE_PARTICLES
 #include "ERF_ParticleData.H"
@@ -1254,6 +1255,8 @@ private:
         amrex::ParallelDescriptor::Barrier("ERF::setRecordSampleLineInfo");
     }
 
+    std::unique_ptr<SampleData> data_sampler = nullptr;
+
     amrex::Vector<std::unique_ptr<std::fstream> > datalog;
     amrex::Vector<std::string> datalogname;
 

Source/ERF.cpp

@@ -1138,6 +1138,12 @@ ERF::InitData_post ()
 
     }
 
+    // DEBUG
+    data_sampler = std::make_unique<SampleData>();
+    data_sampler->get_sample_data(vars_new);
+    data_sampler->write_sample_data(geom);
+    exit(0);
+
     if (pp.contains("sample_line_log") && pp.contains("sample_line"))
     {
         int lev = 0;

Source/IO/ERF_SampleData.H

@@ -0,0 +1,184 @@
+#ifndef ERF_SAMPLEDATA_H
+#define ERF_SAMPLEDATA_H
+
+#include <memory>
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_MultiFabUtil.H>
+#include <AMReX_PlotFileUtil.H>
+
+#include <ERF_IndexDefines.H>
+
+struct LineSampler
+{
+
+    LineSampler ()
+    {
+        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);
+
+            // Allocate space for level indicator
+            m_lev.resize(n_line_dir,0);
+
+            // Allocate space for MF pointers
+            m_ls_mf.resize(n_line_lo);
+        }
+    }
+
+    void
+    get_line_mfs (amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new)
+    {
+        int nlev  = vars_new.size();
+        int nline = m_bnd_bx.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) {
+                amrex::MultiFab& mf = vars_new[ilev][Vars::cons];
+                  m_lev[iline] = ilev;
+                m_ls_mf[iline] = get_line_data(mf, 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
+    }
+
+    void
+    write_line_mfs (amrex::Vector<amrex::Geometry>& geom)
+    {
+        int lev = m_lev[0];
+
+        amrex::Real time = 0.0;
+        amrex::Vector<int> level_steps  = {0};
+        amrex::Vector<amrex::IntVect> ref_ratio = {amrex::IntVect(1)};
+        amrex::Vector<std::string> varnames = {"r", "rT", "rKE", "rQKE", "rS"};
+        std::string name_base = "plt_line_";
+
+        std::string plotfilename = amrex::Concatenate(name_base, 0, 5);
+
+        amrex::Vector<const amrex::MultiFab*> mfp = {&(m_ls_mf[0])};
+        amrex::Vector<amrex::Geometry> temp = {geom[lev]};
+
+        // HACK GEOM TO 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(3,0);
+        amrex::Box m_dom = m_bnd_bx[0];
+        amrex::RealBox m_rb;
+        for (int i(0); i<3; ++i) {
+            int offset = (i==2) ? 0 : 1;
+            amrex::Real xlo = plo[i] + ( m_dom.smallEnd(i) ) * dx[i];
+            amrex::Real xhi = plo[i] + ( m_dom.bigEnd(i) + offset ) * dx[i];
+
+            m_rb.setLo(i,xlo);
+            m_rb.setHi(i,xhi);
+
+            is_per[i] = geom[lev].isPeriodic(i);
+        }
+
+
+        m_geom[0].define(m_dom, &m_rb, geom[lev].Coord(), is_per.data());
+        WriteMultiLevelPlotfile(plotfilename, 1, mfp,
+                                varnames, m_geom, time, level_steps, ref_ratio);
+
+
+        /*
+        WriteMultiLevelPlotfile(plotfilename, 1, mfp,
+                                varnames, temp, time, level_steps, ref_ratio);
+        */
+    }
+
+    amrex::Vector<int> m_dir;
+    amrex::Vector<int> m_lev;
+    amrex::Vector<amrex::Box> m_bnd_bx;
+    amrex::Vector<amrex::MultiFab> m_ls_mf;
+};
+
+class SampleData
+{
+public:
+    explicit SampleData ()
+    {
+        amrex::ParmParse pp("erf");
+        if(pp.contains("sample_line")) m_ls = std::make_unique<LineSampler>();
+    }
+
+    void
+    get_sample_data (amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new)
+    {
+        if (m_ls) m_ls->get_line_mfs(vars_new);
+    }
+
+    void
+    write_sample_data (amrex::Vector<amrex::Geometry>& geom)
+    {
+        if (m_ls) m_ls->write_line_mfs(geom);
+    }
+
+private:
+
+    // Geometry objects for all levels
+    amrex::Vector<amrex::Geometry> m_geom;
+
+    // Variables for IO
+    amrex::Vector<amrex::Vector<std::string>> m_var_names;
+
+    // Structures for getting line MFs
+    std::unique_ptr<LineSampler> m_ls = nullptr;
+
+    // Structures for getting plane MFs
+    //std::unique_ptr<PlaneSampler> m_ps = nullptr;
+};
+#endif

Source/IO/ERF_WriteScalarProfiles.cpp

@@ -302,8 +302,14 @@ ERF::sample_lines (int lev, Real time, IntVect cell, MultiFab& mf)
     // In this case we sample in the vertical direction so dir = 2
     // The "k" value of "cell" is ignored
     //
+    AllPrint() << "PRE GET LINE\n";
+    IntVect cell2 = cell; cell2[2] += 20;
+    IntVect cell3 = cell; cell3[0] += 20;
+    Box bnd_bx_z(cell,cell2);
+    Box bnd_bx_x(cell,cell3);
     int dir = 2;
-    MultiFab my_line       = get_line_data(mf,              dir, cell);
+    MultiFab my_line       = get_line_data(mf,              dir, cell, bnd_bx_z);
+    MultiFab my_line2      = get_line_data(mf,              0  , cell, bnd_bx_x);
     MultiFab my_line_vels  = get_line_data(mf_vels,         dir, cell);
     MultiFab my_line_tau11 = get_line_data(*Tau11_lev[lev], dir, cell);
     MultiFab my_line_tau12 = get_line_data(*Tau12_lev[lev], dir, cell);
@@ -312,6 +318,19 @@ ERF::sample_lines (int lev, Real time, IntVect cell, MultiFab& mf)
     MultiFab my_line_tau23 = get_line_data(*Tau23_lev[lev], dir, cell);
     MultiFab my_line_tau33 = get_line_data(*Tau33_lev[lev], dir, cell);
 
+    Vector<Real> low  = {  16.0,  16.0, 80.0};
+    Vector<Real> high = {1008.0, 512.0, 80.0};
+    RealBox bnd_rbx(low.data(),high.data());
+    std::unique_ptr<MultiFab> test = get_slice_data(2, 80.0, mf, geom[lev], 0, 1, true, bnd_rbx);
+
+    AllPrint() << "CHECK Z: " << bnd_bx_z << ' ' << my_line.boxArray() << ' '
+               << my_line.DistributionMap() << "\n";
+    AllPrint() << "CHECK X: " << bnd_bx_x << ' ' << my_line2.boxArray() << ' '
+               << my_line2.DistributionMap() << "\n";
+    AllPrint() << "CHECK T: " << bnd_rbx << ' ' << test->boxArray() << ' '
+               << test->DistributionMap() << "\n";
+    exit(0);
+
     for (MFIter mfi(my_line, false); mfi.isValid(); ++mfi)
     {
         // HERE DO WHATEVER YOU WANT TO THE DATA BEFORE WRITING

Source/IO/Make.package

@@ -14,6 +14,8 @@ CEXE_sources += ERF_WriteScalarProfiles.cpp
 
 CEXE_sources += ERF_console_io.cpp
 
+CEXE_headers += ERF_SampleData.H
+
 ifeq ($(USE_NETCDF), TRUE)
   CEXE_sources += ERF_ReadFromWRFBdy.cpp
   CEXE_sources += ERF_ReadFromWRFInput.cpp