From 4a91a18cd9ff75a04fc3c6fa47b9b1f88fb255f4 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Thu, 17 Oct 2024 18:16:49 -0700 Subject: [PATCH 1/5] Super rough but it plotted a line. --- Exec/ABL/inputs_sample | 8 +- Source/ERF.H | 3 + Source/ERF.cpp | 6 + Source/IO/ERF_SampleData.H | 184 ++++++++++++++++++++++++++ Source/IO/ERF_WriteScalarProfiles.cpp | 21 ++- Source/IO/Make.package | 2 + 6 files changed, 221 insertions(+), 3 deletions(-) create mode 100644 Source/IO/ERF_SampleData.H diff --git a/Exec/ABL/inputs_sample b/Exec/ABL/inputs_sample index bc5659681..5c9f52852 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,7 +28,11 @@ 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 diff --git a/Source/ERF.H b/Source/ERF.H index e9ed7ad52..0b1085bf3 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -40,6 +40,7 @@ #include #include #include +#include #ifdef ERF_USE_PARTICLES #include "ERF_ParticleData.H" @@ -1254,6 +1255,8 @@ private: amrex::ParallelDescriptor::Barrier("ERF::setRecordSampleLineInfo"); } + std::unique_ptr data_sampler = nullptr; + amrex::Vector > datalog; amrex::Vector datalogname; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 427582722..32897e9b7 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1138,6 +1138,12 @@ 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; diff --git a/Source/IO/ERF_SampleData.H b/Source/IO/ERF_SampleData.H new file mode 100644 index 000000000..88bcfea01 --- /dev/null +++ b/Source/IO/ERF_SampleData.H @@ -0,0 +1,184 @@ +#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(); + + // Loop over each line + for (int iline(0); iline=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& 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_"; + + std::string plotfilename = amrex::Concatenate(name_base, 0, 5); + + amrex::Vector mfp = {&(m_ls_mf[0])}; + amrex::Vector temp = {geom[lev]}; + + // HACK GEOM TO LINE + auto plo = geom[lev].ProbLo(); + auto dx = geom[lev].CellSize(); + + 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]; + + 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 m_dir; + amrex::Vector m_lev; + amrex::Vector m_bnd_bx; + amrex::Vector m_ls_mf; +}; + +class SampleData +{ +public: + explicit SampleData () + { + amrex::ParmParse pp("erf"); + if(pp.contains("sample_line")) m_ls = std::make_unique(); + } + + void + get_sample_data (amrex::Vector>& vars_new) + { + if (m_ls) m_ls->get_line_mfs(vars_new); + } + + void + write_sample_data (amrex::Vector& geom) + { + if (m_ls) m_ls->write_line_mfs(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; +}; +#endif diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index b531afc85..261f879ea 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/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 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 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 From 95f5afd950fbcd067acb3a2a55d2bf165c288dab Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Fri, 18 Oct 2024 13:40:59 -0700 Subject: [PATCH 2/5] Need to test on multi-proc. --- Exec/ABL/inputs_DataSampler | 74 +++++++ Exec/ABL/inputs_sample | 8 +- Source/ERF.H | 3 + Source/ERF.cpp | 21 +- Source/IO/ERF_SampleData.H | 265 +++++++++++++++++++++----- Source/IO/ERF_WriteScalarProfiles.cpp | 21 +- 6 files changed, 310 insertions(+), 82 deletions(-) create mode 100644 Exec/ABL/inputs_DataSampler 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 From db97e5ad7d8142fa580771e7fbf3b501307c7484 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Fri, 18 Oct 2024 13:52:30 -0700 Subject: [PATCH 3/5] codespell fix. --- Source/ERF.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index b13ccef20..6582f146d 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1173,7 +1173,7 @@ ERF::InitData_post () } - // Create object to do line and plane sampling if neeeded + // 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); } From 366e5bbbd79f590b242194eba7f80d9e533367fa Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Mon, 21 Oct 2024 14:03:59 -0700 Subject: [PATCH 4/5] Plot T and Wsp and clean up the realBox location. --- Docs/sphinx_doc/Inputs.rst | 68 +++++++++++++++++++++++++++++ Source/IO/ERF_SampleData.H | 89 ++++++++++++++++++++++++++++++++------ 2 files changed, 144 insertions(+), 13 deletions(-) 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/Source/IO/ERF_SampleData.H b/Source/IO/ERF_SampleData.H index 6a1137a83..30ddd9454 100644 --- a/Source/IO/ERF_SampleData.H +++ b/Source/IO/ERF_SampleData.H @@ -71,6 +71,7 @@ struct LineSampler { 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) { - amrex::MultiFab& mf = vars_new[ilev][Vars::cons]; + + // 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, dir, cell, bnd_bx); + 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(); @@ -99,7 +131,7 @@ struct LineSampler amrex::Vector& geom) { std::string name_base = "plt_line_"; - amrex::Vector varnames = {"r", "rT", "rKE", "rQKE", "rS"}; + amrex::Vector varnames = {"T", "Wsp"}; int nline = m_ls_mf.size(); for (int iline(0); iline=0; --ilev) { - amrex::MultiFab& mf = vars_new[ilev][Vars::cons]; + + // 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, geom[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 @@ -262,7 +325,7 @@ struct PlaneSampler amrex::Vector& geom) { std::string name_base = "plt_plane_"; - amrex::Vector varnames = {"r", "rT", "rKE", "rQKE", "rS"}; + amrex::Vector varnames = {"T", "Wsp"}; int nplane = m_ps_mf.size(); for (int iplane(0); iplane is_per(AMREX_SPACEDIM,0); for (int d(0); d>& vars_new) { if (m_ls) m_ls->get_line_mfs(vars_new); - if (m_ps) m_ps->get_plane_mfs(geom,vars_new); + if (m_ps) m_ps->get_plane_mfs(geom, vars_new); } void From 808976c49a574ee749b9462f5d7f50fb061ccf5d Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Mon, 21 Oct 2024 14:11:53 -0700 Subject: [PATCH 5/5] Bump AMReX hash. --- Submodules/AMReX | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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