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/Exec/ABL/inputs_sample b/Exec/ABL/inputs_sample index 5c9f52852..bc5659681 100644 --- a/Exec/ABL/inputs_sample +++ b/Exec/ABL/inputs_sample @@ -20,7 +20,7 @@ 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 @@ -28,11 +28,7 @@ erf.sample_point_log = my_sample_point erf.sample_point = 6 63 5 erf.sample_line_log = my_sample_line -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 +erf.sample_line = 63 63 5 # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed diff --git a/Source/ERF.H b/Source/ERF.H index 0b1085bf3..999ba0be4 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -1255,6 +1255,9 @@ 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; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 32897e9b7..90986eef6 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -539,6 +539,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) ) { @@ -1138,12 +1144,6 @@ ERF::InitData_post () } - // DEBUG - data_sampler = std::make_unique(); - 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; @@ -1174,6 +1174,11 @@ ERF::InitData_post () } + // Create object to do line and plane sampling if neeeded + 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); @@ -1501,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 index 88bcfea01..6a1137a83 100644 --- a/Source/IO/ERF_SampleData.H +++ b/Source/IO/ERF_SampleData.H @@ -12,7 +12,6 @@ struct LineSampler { - LineSampler () { amrex::ParmParse pp("erf"); @@ -30,7 +29,7 @@ struct LineSampler 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++) { + 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]); @@ -41,7 +40,7 @@ struct LineSampler 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++) { + 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]); @@ -87,98 +86,264 @@ struct LineSampler // 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; + if (bnd_bx == min_bnd_bx) { continue; } } // ilev }// iline } void - write_line_mfs (amrex::Vector& geom) + write_line_mfs (amrex::Vector& time, + amrex::Vector& level_steps, + amrex::Vector& ref_ratio, + amrex::Vector& geom) { - int lev = m_lev[0]; - - amrex::Real time = 0.0; - amrex::Vector level_steps = {0}; - amrex::Vector ref_ratio = {amrex::IntVect(1)}; - amrex::Vector varnames = {"r", "rT", "rKE", "rQKE", "rS"}; std::string name_base = "plt_line_"; + amrex::Vector varnames = {"r", "rT", "rKE", "rQKE", "rS"}; + + 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; +}; - std::string plotfilename = amrex::Concatenate(name_base, 0, 5); - amrex::Vector mfp = {&(m_ls_mf[0])}; - amrex::Vector temp = {geom[lev]}; +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); + } - // HACK GEOM TO LINE - auto plo = geom[lev].ProbLo(); - auto dx = geom[lev].CellSize(); + // 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; + } - amrex::Vector m_geom; m_geom.resize(1); - amrex::Vector 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]; + // Parse directionality + m_dir.resize(n_plane_dir); + pp.queryarr("sample_plane_dir",m_dir,0,n_plane_dir); - m_rb.setLo(i,xlo); - m_rb.setHi(i,xhi); + // Allocate space for level indicator + m_lev.resize(n_plane_dir,0); - is_per[i] = geom[lev].isPeriodic(i); + // 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 = vars_new[0][Vars::cons].nComp(); + bool interpolate = true; + + // Loop over each plane + for (int iplane(0); iplane=0; --ilev) { + amrex::MultiFab& mf = vars_new[ilev][Vars::cons]; + m_lev[iplane] = ilev; + m_ps_mf[iplane] = get_slice_data(dir, point, mf, 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; } - 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); + } // 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 = {"r", "rT", "rKE", "rQKE", "rS"}; - /* - WriteMultiLevelPlotfile(plotfilename, 1, mfp, - varnames, temp, time, level_steps, ref_ratio); - */ + 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_bx; - amrex::Vector m_ls_mf; + amrex::Vector m_bnd_rbx; + amrex::Vector> m_ps_mf; }; + class SampleData { public: - explicit SampleData () + explicit SampleData (bool do_line=false, + bool do_plane=false) { - amrex::ParmParse pp("erf"); - if(pp.contains("sample_line")) m_ls = std::make_unique(); + if(do_line) m_ls = std::make_unique(); + if(do_plane) m_ps = std::make_unique(); } void - get_sample_data (amrex::Vector>& vars_new) + 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& geom) + 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(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: - // Geometry objects for all levels - amrex::Vector m_geom; - - // Variables for IO - amrex::Vector> m_var_names; - // Structures for getting line MFs std::unique_ptr m_ls = nullptr; // Structures for getting plane MFs - //std::unique_ptr m_ps = nullptr; + std::unique_ptr m_ps = nullptr; }; #endif diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index 261f879ea..b531afc85 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/Source/IO/ERF_WriteScalarProfiles.cpp @@ -302,14 +302,8 @@ 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, bnd_bx_z); - MultiFab my_line2 = get_line_data(mf, 0 , cell, bnd_bx_x); + MultiFab my_line = get_line_data(mf, dir, cell); 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); @@ -318,19 +312,6 @@ 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 low = { 16.0, 16.0, 80.0}; - Vector high = {1008.0, 512.0, 80.0}; - RealBox bnd_rbx(low.data(),high.data()); - std::unique_ptr 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