diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 791e4d1b9..b3559ff68 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -718,6 +718,74 @@ The requested output files have the following columns: #. SGS turbulence dissipation, :math:`\epsilon` (m2/s3) +Data Sampling Outputs +================== + +Data along query lines or planes may be output during the simulation if +``erf.do_line_sampling = true`` or ``erf.do_plane_sampling = true``, respectively. +The potential temperature and wind-speed will be written to native ``plt_line/plane`` +at the step frequency dictated by ``erf.sampler_interval = ``. For line sampling, +users must prescribe ``sample_line_lo`` and ``sample_line_hi`` inputs which are 3 integer +values corresponding to the (i,j,k) indices at the beginning and end of the line. +Additionally, users must specify ``sample_line_dir`` to prescribed the direction of +the line. The same inputs are used for the plane sampling except that ``sample_plane_lo/hi`` +must be the physical locations of the plane corners. This output functionality has +not been implemented for terrain. + +.. _list-of-parameters-10b: + + +List of Parameters +------------------ + ++-------------------------------+------------------+----------------+----------------+ +| Parameter | Definition | Acceptable | Default | +| | | Values | | ++===============================+==================+================+================+ +| **erf.sampler_interval** | Output | Integer | -1 | +| | frequency | | | ++-------------------------------+------------------+----------------+----------------+ +| **erf.do_line_sampling** | Flag to do line | Boolean | false | +| | sampling | | | +| | | | | ++-------------------------------+------------------+----------------+----------------+ +| **erf.do_plane_sampling** | Flag to do plane | Boolean | false | +| | sampling | | | +| | | | | ++-------------------------------+------------------+----------------+----------------+ +| **erf.sample_line_dir** | Directionality | Integer | None | +| | of the line | | | ++-------------------------------+------------------+----------------+----------------+ +| **erf.sample_plane_dir** | Directionality | Integer | None | +| | of the plane | | | ++-------------------------------+------------------+----------------+----------------+ +| **erf.sample_line_lo/hi** | Bounding (i,j,k) | 3 Integers per | None | +| | on the line(s) | line | | ++-------------------------------+------------------+----------------+----------------+ +| **erf.sample_plane_lo/hi** | Bounding point | 3 Reals per | None | +| | on the plane(s) | plane | | ++-------------------------------+------------------+----------------+----------------+ + +.. _examples-of-usage-10b: + +Example of Usage +----------------- + +:: + + erf.sampler_interval = 1 # Write plt files every step + + erf.do_line_sampling = true # Do line sampling + erf.sample_line_lo = 5 32 5 10 32 5 # Lo points for two lines + erf.sample_line_hi = 5 32 25 1000 32 5 # Hi points for two lines + erf.sample_line_dir = 2 0 # One line in z and one in x + + erf.do_plane_sampling = true # Do plane sampling + erf.sample_plane_lo = 48.0 48.0 32.0 # Lo points for one plane + erf.sample_plane_hi = 320.0 320.0 32.0 # Hi points for one plane + erf.sample_plane_dir = 2 # One plane with z normal + + Advection Schemes ================= diff --git a/Exec/ABL/inputs_DataSampler b/Exec/ABL/inputs_DataSampler new file mode 100644 index 000000000..ac46484c4 --- /dev/null +++ b/Exec/ABL/inputs_DataSampler @@ -0,0 +1,74 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 4000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 1024 1024 1024 +amr.n_cell = 64 64 64 + +geometry.is_periodic = 1 1 0 + +zlo.type = "NoSlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +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 + +erf.data_log = my_data_file + +erf.sampler_interval = 1 +erf.do_line_sampling = true +erf.sample_line_lo = 5 32 5 10 32 5 +erf.sample_line_hi = 5 32 25 1000 32 5 +erf.sample_line_dir = 2 0 + +erf.do_plane_sampling = true +erf.sample_plane_lo = 48.0 48.0 32.0 +erf.sample_plane_hi = 320.0 320.0 32.0 +erf.sample_plane_dir = 2 + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 10 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 + +erf.init_type = "uniform" + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + +# Higher values of perturbations lead to instability +# Instability seems to be coming from BC +prob.U_0_Pert_Mag = 0.08 +prob.V_0_Pert_Mag = 0.08 # +prob.W_0_Pert_Mag = 0.0 diff --git a/Source/ERF.H b/Source/ERF.H index 26c1f01a4..a977fa35a 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -44,6 +44,7 @@ #include #include #include +#include #ifdef ERF_USE_PARTICLES #include "ERF_ParticleData.H" @@ -1263,6 +1264,11 @@ private: amrex::ParallelDescriptor::Barrier("ERF::setRecordSampleLineInfo"); } + // Data sampler for line and plane output + int sampler_interval = -1; + amrex::Real sampler_per = -1.0; + std::unique_ptr data_sampler = nullptr; + amrex::Vector > datalog; amrex::Vector datalogname; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 85ffb5aec..ef6ec5bfc 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -538,6 +538,12 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) } } + // Write plane/line sampler data + if (is_it_time_for_action(nstep, time, dt_lev0, sampler_interval, sampler_per) && (data_sampler) ) { + data_sampler->get_sample_data(geom, vars_new); + data_sampler->write_sample_data(t_new, istep, ref_ratio, geom); + } + // Moving terrain if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) { @@ -1167,6 +1173,11 @@ ERF::InitData_post () } + // Create object to do line and plane sampling if needed + bool do_line = false; bool do_plane = false; + pp.query("do_line_sampling",do_line); pp.query("do_plane_sampling",do_plane); + if (do_line || do_plane) { data_sampler = std::make_unique(do_line, do_plane); } + #ifdef ERF_USE_EB bool write_eb_surface = false; pp.query("write_eb_surface", write_eb_surface); @@ -1495,6 +1506,10 @@ ERF::ReadParameters () pp.query("column_loc_y", column_loc_y); pp.query("column_file_name", column_file_name); + // Sampler output frequency + pp.query("sampler_per", sampler_per); + pp.query("sampler_interval", sampler_interval); + // Specify information about outputting planes of data pp.query("output_bndry_planes", output_bndry_planes); pp.query("bndry_output_planes_interval", bndry_output_planes_interval); diff --git a/Source/IO/ERF_SampleData.H b/Source/IO/ERF_SampleData.H new file mode 100644 index 000000000..30ddd9454 --- /dev/null +++ b/Source/IO/ERF_SampleData.H @@ -0,0 +1,412 @@ +#ifndef ERF_SAMPLEDATA_H +#define ERF_SAMPLEDATA_H + +#include + +#include +#include +#include +#include + +#include + +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 idx_lo; idx_lo.resize(n_line_lo*AMREX_SPACEDIM); + amrex::Vector 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 idx_hi; idx_hi.resize(n_line_hi*AMREX_SPACEDIM); + amrex::Vector 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>& vars_new) + { + int nlev = vars_new.size(); + int nline = m_bnd_bx.size(); + int ncomp = 2; + + // Loop over each line + for (int iline(0); iline=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{&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); +#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& tfab = vars_new[ilev][Vars::cons].array(mfi); + auto const& wfab = mf_cc_vel.array(mfi); + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + dfab(i,j,k,0) = tfab(i,j,k,1)/tfab(i,j,k,0); + dfab(i,j,k,1) = std::sqrt(wfab(i,j,k,0)*wfab(i,j,k,0) + + wfab(i,j,k,1)*wfab(i,j,k,1) + + wfab(i,j,k,2)*wfab(i,j,k,2)) ; + }); + + } + + 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 + } + + void + write_line_mfs (amrex::Vector& time, + amrex::Vector& level_steps, + amrex::Vector& ref_ratio, + amrex::Vector& geom) + { + std::string name_base = "plt_line_"; + amrex::Vector varnames = {"T", "Wsp"}; + + int nline = m_ls_mf.size(); + for (int iline(0); iline m_level_steps = {level_steps[lev]}; + amrex::Vector 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 m_geom; m_geom.resize(1); + amrex::Vector is_per(AMREX_SPACEDIM,0); + amrex::Box m_dom = m_bnd_bx[iline]; + amrex::RealBox m_rb; + for (int d(0); d mf = {&(m_ls_mf[iline])}; + + // Write each line + WriteMultiLevelPlotfile(plotfilename, 1, mf, + varnames, m_geom, m_time, + m_level_steps, m_ref_ratio); + } + } + + amrex::Vector m_dir; + amrex::Vector m_lev; + amrex::Vector m_bnd_bx; + amrex::Vector m_ls_mf; +}; + + +struct PlaneSampler +{ + PlaneSampler () + { + 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 r_lo; r_lo.resize(n_plane_lo*AMREX_SPACEDIM); + amrex::Vector> 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 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 r_hi; r_hi.resize(n_plane_hi*AMREX_SPACEDIM); + amrex::Vector> 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 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); + + // Allocate space for level indicator + m_lev.resize(n_plane_dir,0); + + // Allocate space for MF pointers + m_ps_mf.resize(n_plane_lo); + } + } + + // This must match what is in AMReX_MultiFabUtil.H + amrex::Box + getIndexBox (const amrex::RealBox& real_box, + const amrex::Geometry& geom) { + amrex::IntVect slice_lo, slice_hi; + + AMREX_D_TERM(slice_lo[0]=static_cast(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));, + slice_lo[1]=static_cast(std::floor((real_box.lo(1) - geom.ProbLo(1))/geom.CellSize(1)));, + slice_lo[2]=static_cast(std::floor((real_box.lo(2) - geom.ProbLo(2))/geom.CellSize(2)));); + + AMREX_D_TERM(slice_hi[0]=static_cast(std::floor((real_box.hi(0) - geom.ProbLo(0))/geom.CellSize(0)));, + slice_hi[1]=static_cast(std::floor((real_box.hi(1) - geom.ProbLo(1))/geom.CellSize(1)));, + slice_hi[2]=static_cast(std::floor((real_box.hi(2) - geom.ProbLo(2))/geom.CellSize(2)));); + + return amrex::Box(slice_lo, slice_hi) & geom.Domain(); + } + + void + get_plane_mfs (amrex::Vector& geom, + amrex::Vector>& vars_new) + { + int nlev = vars_new.size(); + int nplane = m_bnd_rbx.size(); + int ncomp = 2; + bool interpolate = true; + + // Loop over each plane + for (int iplane(0); iplane=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{&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); +#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& tfab = vars_new[ilev][Vars::cons].array(mfi); + auto const& wfab = mf_cc_vel.array(mfi); + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + dfab(i,j,k,0) = tfab(i,j,k,1)/tfab(i,j,k,0); + dfab(i,j,k,1) = std::sqrt(wfab(i,j,k,0)*wfab(i,j,k,0) + + wfab(i,j,k,1)*wfab(i,j,k,1) + + wfab(i,j,k,2)*wfab(i,j,k,2)) ; + }); + + } + + 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) { continue; } + + } // ilev + }// iplane + } + + void + write_plane_mfs (amrex::Vector& time, + amrex::Vector& level_steps, + amrex::Vector& ref_ratio, + amrex::Vector& geom) + { + std::string name_base = "plt_plane_"; + amrex::Vector varnames = {"T", "Wsp"}; + + int nplane = m_ps_mf.size(); + for (int iplane(0); iplane m_level_steps = {level_steps[lev]}; + amrex::Vector 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 is_per(AMREX_SPACEDIM,0); + for (int d(0); d 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 = amrex::Concatenate(name_base, iplane , 5); + name_plane += "_step_"; + std::string plotfilename = amrex::Concatenate(name_plane, m_level_steps[0], 5); + + // Get the data + amrex::Vector mf = {m_ps_mf[iplane].get()}; + + // Write each plane + WriteMultiLevelPlotfile(plotfilename, 1, mf, + varnames, m_geom, m_time, + m_level_steps, m_ref_ratio); + } // iplane + } + + amrex::Vector m_dir; + amrex::Vector m_lev; + amrex::Vector m_bnd_rbx; + amrex::Vector> m_ps_mf; +}; + + +class SampleData +{ +public: + explicit SampleData (bool do_line=false, + bool do_plane=false) + { + if(do_line) m_ls = std::make_unique(); + if(do_plane) m_ps = std::make_unique(); + } + + void + get_sample_data (amrex::Vector& geom, + amrex::Vector>& vars_new) + { + if (m_ls) m_ls->get_line_mfs(vars_new); + if (m_ps) m_ps->get_plane_mfs(geom, vars_new); + } + + void + write_sample_data (amrex::Vector& time, + amrex::Vector& level_steps, + amrex::Vector& ref_ratio, + amrex::Vector& geom) + { + if (m_ls) m_ls->write_line_mfs(time, level_steps, ref_ratio, geom); + if (m_ls) m_ps->write_plane_mfs(time, level_steps, ref_ratio, geom); + } + +private: + + // Structures for getting line MFs + std::unique_ptr m_ls = nullptr; + + // Structures for getting plane MFs + std::unique_ptr m_ps = nullptr; +}; +#endif diff --git a/Source/IO/Make.package b/Source/IO/Make.package index 4f6f079c9..d2369a85f 100644 --- a/Source/IO/Make.package +++ b/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 diff --git a/Submodules/AMReX b/Submodules/AMReX index 62c2a81ea..c333708d5 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit 62c2a81eac7862d526e5861ef2befc00b7f5b759 +Subproject commit c333708d59f01ab664363ba426e85ad24c1fb23d