ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_SampleData.H
Go to the documentation of this file.
1 /**
2  * \file ERF_SampleData.H
3  */
4 #ifndef ERF_SAMPLEDATA_H
5 #define ERF_SAMPLEDATA_H
6 
7 #include <memory>
8 
9 #include <AMReX_ParmParse.H>
10 #include <AMReX_MultiFab.H>
11 #include <AMReX_MultiFabUtil.H>
12 #include <AMReX_PlotFileUtil.H>
13 
14 #include <ERF_IndexDefines.H>
15 #include <ERF_EOS.H>
16 #include <ERF_Container.H>
17 
18 #include <iostream>
19 
20 struct LineSampler
21 {
22  LineSampler ()
23  {
24  amrex::ParmParse pp("erf");
25 
26  bool has_line_idx = pp.contains("sample_line_lo") || pp.contains("sample_line_hi");
27  bool has_line_real = pp.contains("sample_line_lo_real") || pp.contains("sample_line_hi_real");
28  if (has_line_idx && has_line_real) {
29  amrex::Abort("Specify only one of erf.sample_line_lo/hi or erf.sample_line_lo_real/hi_real");
30  }
31 
32  // Count number of lo and hi points define the line
33  int n_line_lo = pp.countval("sample_line_lo") / AMREX_SPACEDIM;
34  int n_line_hi = pp.countval("sample_line_hi") / AMREX_SPACEDIM;
35  int n_line_lo_real = pp.countval("sample_line_lo_real") / AMREX_SPACEDIM;
36  int n_line_hi_real = pp.countval("sample_line_hi_real") / AMREX_SPACEDIM;
37  int n_line_dir = pp.countval("sample_line_dir");
38  if (has_line_idx) {
39  AMREX_ALWAYS_ASSERT( (n_line_lo==n_line_hi ) &&
40  (n_line_lo==n_line_dir) );
41  } else if (has_line_real) {
42  AMREX_ALWAYS_ASSERT( (n_line_lo_real==n_line_hi_real ) &&
43  (n_line_lo_real==n_line_dir) );
44  }
45 
46  // Parse the data
47  int nline = 0;
48  if (has_line_idx) {
49  nline = n_line_lo;
50  if (nline > 0) {
51  // Parse lo
52  amrex::Vector<int> idx_lo; idx_lo.resize(nline*AMREX_SPACEDIM);
53  amrex::Vector<amrex::IntVect> iv_lo; iv_lo.resize(nline);
54  pp.queryarr("sample_line_lo",idx_lo,0,nline*AMREX_SPACEDIM);
55  for (int i(0); i < nline; i++) {
56  amrex::IntVect iv(idx_lo[AMREX_SPACEDIM*i+0],
57  idx_lo[AMREX_SPACEDIM*i+1],
58  idx_lo[AMREX_SPACEDIM*i+2]);
59  iv_lo[i] = iv;
60  }
61 
62  // Parse hi
63  amrex::Vector<int> idx_hi; idx_hi.resize(nline*AMREX_SPACEDIM);
64  amrex::Vector<amrex::IntVect> iv_hi; iv_hi.resize(nline);
65  pp.queryarr("sample_line_hi",idx_hi,0,nline*AMREX_SPACEDIM);
66  for (int i(0); i < nline; i++) {
67  amrex::IntVect iv(idx_hi[AMREX_SPACEDIM*i+0],
68  idx_hi[AMREX_SPACEDIM*i+1],
69  idx_hi[AMREX_SPACEDIM*i+2]);
70  iv_hi[i] = iv;
71  }
72 
73  // Construct vector of bounding boxes
74  m_bnd_bx.resize(nline);
75  for (int i = 0; i < nline; i++){
76  amrex::Box lbx(iv_lo[i],iv_hi[i]);
77  m_bnd_bx[i] = lbx;
78  }
79  }
80  } else if (has_line_real) {
81  m_use_real_bx = true;
82  nline = n_line_lo_real;
83  if (nline > 0) {
84  // Parse lo
85  amrex::Vector<amrex::Real> real_lo; real_lo.resize(nline*AMREX_SPACEDIM);
86  amrex::Vector<amrex::Vector<amrex::Real>> rv_lo;
87  pp.queryarr("sample_line_lo_real",real_lo,0,nline*AMREX_SPACEDIM);
88  for (int i(0); i < nline; i++) {
89  amrex::Vector<amrex::Real> rv = {real_lo[AMREX_SPACEDIM*i+0],
90  real_lo[AMREX_SPACEDIM*i+1],
91  real_lo[AMREX_SPACEDIM*i+2]};
92  rv_lo.push_back(rv);
93  }
94 
95  // Parse hi
96  amrex::Vector<amrex::Real> real_hi; real_hi.resize(nline*AMREX_SPACEDIM);
97  amrex::Vector<amrex::Vector<amrex::Real>> rv_hi;
98  pp.queryarr("sample_line_hi_real",real_hi,0,nline*AMREX_SPACEDIM);
99  for (int i(0); i < nline; i++) {
100  amrex::Vector<amrex::Real> rv = {real_hi[AMREX_SPACEDIM*i+0],
101  real_hi[AMREX_SPACEDIM*i+1],
102  real_hi[AMREX_SPACEDIM*i+2]};
103  rv_hi.push_back(rv);
104  }
105 
106  // Construct vector of bounding real boxes
107  m_bnd_rbx.resize(nline);
108  for (int i = 0; i < nline; i++){
109  amrex::RealBox rbx(rv_lo[i].data(),rv_hi[i].data());
110  m_bnd_rbx[i] = rbx;
111  }
112 
113  // Construct vector of bounding boxes
114  m_bnd_bx.resize(nline);
115  }
116  }
117 
118  if (nline > 0) {
119  // Parse directionality
120  m_dir.resize(n_line_dir);
121  pp.queryarr("sample_line_dir",m_dir,0,n_line_dir);
122 
123  // Parse names
124  std::string name_base = "plt_line_";
125  m_name.resize(nline);
126  int n_names = pp.countval("sample_line_name");
127  if (n_names > 0) {
128  AMREX_ALWAYS_ASSERT( n_names==nline );
129  pp.queryarr("sample_line_name",m_name,0,n_names);
130  } else {
131  for (int iline(0); iline<nline; ++iline) {
132  m_name[iline] = amrex::Concatenate(name_base, iline , 5);
133  }
134  }
135 
136  // Allocate space for level indicator
137  m_lev.resize(n_line_dir,0);
138 
139  // Allocate space for MF pointers
140  m_ls_mf.resize(nline);
141 
142  // Get requested vars
143  if (pp.countval("line_sampling_vars") > 0) {
144  m_varnames.clear();
145  amrex::Vector<std::string> requested_vars;
146  pp.queryarr("line_sampling_vars",requested_vars);
147  amrex::Print() << "Selected line sampling vars :";
148  if (containerHasElement(requested_vars, "density")) {
149  m_varnames.push_back("density");
150  amrex::Print() << " " << "density";
151  }
152  if (containerHasElement(requested_vars, "x_velocity")) {
153  m_varnames.push_back("x_velocity");
154  amrex::Print() << " " << "x_velocity";
155  }
156  if (containerHasElement(requested_vars, "y_velocity")) {
157  m_varnames.push_back("y_velocity");
158  amrex::Print() << " " << "y_velocity";
159  }
160  if (containerHasElement(requested_vars, "z_velocity")) {
161  m_varnames.push_back("z_velocity");
162  amrex::Print() << " " << "z_velocity";
163  }
164  if (containerHasElement(requested_vars, "magvel")) {
165  m_varnames.push_back("magvel");
166  amrex::Print() << " " << "magvel";
167  }
168  if (containerHasElement(requested_vars, "theta")) {
169  m_varnames.push_back("theta");
170  amrex::Print() << " " << "theta";
171  }
172  if (containerHasElement(requested_vars, "qv")) {
173  m_varnames.push_back("qv");
174  amrex::Print() << " " << "qv";
175  }
176  if (containerHasElement(requested_vars, "qc")) {
177  m_varnames.push_back("qc");
178  amrex::Print() << " " << "qc";
179  }
180  if (containerHasElement(requested_vars, "pressure")) {
181  m_varnames.push_back("pressure");
182  amrex::Print() << " " << "pressure";
183  }
184  if (containerHasElement(requested_vars, "sgs_tke")) {
185  m_varnames.push_back("sgs_tke");
186  amrex::Print() << " " << "sgs_tke";
187  }
188  if (containerHasElement(requested_vars, "sgs_tau13")) {
189  m_varnames.push_back("sgs_tau13");
190  amrex::Print() << " " << "sgs_tau13";
191  }
192  if (containerHasElement(requested_vars, "sgs_tau23")) {
193  m_varnames.push_back("sgs_tau23");
194  amrex::Print() << " " << "sgs_tau23";
195  }
196  if (containerHasElement(requested_vars, "sgs_hfx3")) {
197  m_varnames.push_back("sgs_hfx3");
198  amrex::Print() << " " << "sgs_hfx3";
199  }
200  amrex::Print() << std::endl;
201  }
202 
203  // Write outputs to text files, one file per variable, with all
204  // times appended to the same file
205  pp.queryAdd("line_sampling_text_output",m_write_ascii);
206  if (m_write_ascii && amrex::ParallelDescriptor::IOProcessor()) {
207  int nvar = static_cast<int>(m_varnames.size());
208  m_datastream.resize(nline * nvar);
209  int i = 0;
210  for (int iline(0); iline<nline; ++iline) {
211  for (int ivar(0); ivar<nvar; ++ivar) {
212  std::string filename = m_name[iline] + "." + m_varnames[ivar];
213  m_datastream[i] = std::make_unique<std::fstream>();
214  m_datastream[i]->open(filename.c_str(),std::ios::out|std::ios::app);
215  if (!m_datastream[i]->good()) {
216  amrex::FileOpenFailed(filename);
217  }
218  i++;
219  }
220  }
221  }
222  }
223  }
224 
225  amrex::Box
226  getIndexBox (const amrex::RealBox& real_box,
227  const amrex::Geometry& geom) {
228  amrex::IntVect slice_lo, slice_hi;
229 
230  AMREX_D_TERM(slice_lo[0]=static_cast<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
231  slice_lo[1]=static_cast<int>(std::floor((real_box.lo(1) - geom.ProbLo(1))/geom.CellSize(1)));,
232  slice_lo[2]=static_cast<int>(std::floor((real_box.lo(2) - geom.ProbLo(2))/geom.CellSize(2))););
233 
234  AMREX_D_TERM(slice_hi[0]=static_cast<int>(std::floor((real_box.hi(0) - geom.ProbLo(0))/geom.CellSize(0)));,
235  slice_hi[1]=static_cast<int>(std::floor((real_box.hi(1) - geom.ProbLo(1))/geom.CellSize(1)));,
236  slice_hi[2]=static_cast<int>(std::floor((real_box.hi(2) - geom.ProbLo(2))/geom.CellSize(2))););
237 
238  return amrex::Box(slice_lo, slice_hi) & geom.Domain();
239  }
240 
241  void
242  write_coords (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& z_phys_cc,
243  amrex::Vector<amrex::Geometry>& geom)
244  {
245  if (!m_write_ascii) return;
246 
247  amrex::Print() << "Writing out line coordinates to text" << std::endl;
248 
249  for (int lev(0); lev < z_phys_cc.size(); ++lev) {
250  // Write one text file per level
251  std::ofstream outfile;
252  if (amrex::ParallelDescriptor::IOProcessor()) {
253  std::string fname = amrex::Concatenate("plt_line_lev", lev, 1);
254  fname += ".zcc";
255  outfile.open(fname);
256 
257  if (!outfile.is_open()) {
258  amrex::AllPrint() << "Could not open " << fname << std::endl;
259  }
260  }
261 
262  // Loop over each line
263  int nline = static_cast<int>(m_ls_mf.size());
264  for (int iline(0); iline<nline; ++iline) {
265  int dir = m_dir[iline];
266  amrex::Box bnd_bx = (m_use_real_bx)
267  ? getIndexBox(m_bnd_rbx[iline], geom[lev])
268  : m_bnd_bx[iline];
269  amrex::IntVect first_cell = bnd_bx.smallEnd();
270 
271  // Create multifab with "sampled" z_phys values
272  amrex::MultiFab line_coords_mf = get_line_data(
273  *z_phys_cc[lev], dir, first_cell, bnd_bx
274  );
275 
276  // Convert multifab to vector
277  amrex::Gpu::HostVector<amrex::Real> vec = sumToLine(
278  line_coords_mf, 0, 1, bnd_bx, dir
279  );
280 
281  // Append to file
282  if (amrex::ParallelDescriptor::IOProcessor()) {
283  for (const auto& zval : vec) {
284  outfile << " " << zval;
285  }
286  outfile << std::endl;
287  }
288  } // line loop
289 
290  outfile.close();
291  } // level loop
292  }
293 
294  void
295  get_sample_data (amrex::Vector<amrex::Geometry>& geom,
296  amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new,
297  const amrex::Vector<amrex::MultiFab*>& tau13_lev = {},
298  const amrex::Vector<amrex::MultiFab*>& tau23_lev = {},
299  const amrex::Vector<amrex::MultiFab*>& hfx3_lev = {},
300  const amrex::Vector<int>& prognostic_tke_available = {})
301  {
302  int nlev = static_cast<int>(vars_new.size());
303  int nline = static_cast<int>(m_bnd_bx.size());
304  int ncomp = static_cast<int>(m_varnames.size());
305 
306  int qv_comp = -1;
307 
308  // Loop over each line
309  for (int iline(0); iline<nline; ++iline) {
310  int dir = m_dir[iline];
311 
312  // Search each level to get the finest data possible
313  for (int ilev(nlev-1); ilev>=0; --ilev) {
314  amrex::Box bnd_bx = (m_use_real_bx)
315  ? getIndexBox(m_bnd_rbx[iline], geom[ilev])
316  : m_bnd_bx[iline];
317  amrex::IntVect cell = bnd_bx.smallEnd();
318 
319  // Construct CC velocities
320  amrex::MultiFab mf_cc_vel;
321  auto ba = vars_new[ilev][Vars::cons].boxArray();
322  auto dm = vars_new[ilev][Vars::cons].DistributionMap();
323  mf_cc_vel.define(ba, dm, AMREX_SPACEDIM, amrex::IntVect(1,1,1));
324  average_face_to_cellcenter(mf_cc_vel,0,
325  amrex::Array<const amrex::MultiFab*,3>{&vars_new[ilev][Vars::xvel],
326  &vars_new[ilev][Vars::yvel],
327  &vars_new[ilev][Vars::zvel]});
328 
329  // Construct MultiFab holding requested variables
330  amrex::MultiFab mf_cc_data;
331  mf_cc_data.define(ba, dm, ncomp, 1);
332 
333  int mf_comp = 0;
334 
335  if (containerHasElement(m_varnames, "density")) {
336  amrex::MultiFab::Copy(mf_cc_data, vars_new[ilev][Vars::cons], Rho_comp, mf_comp, 1, 0);
337  mf_comp += 1;
338  }
339 
340  if (containerHasElement(m_varnames, "x_velocity")) {
341  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 0, mf_comp, 1, 0);
342  mf_comp += 1;
343  }
344  if (containerHasElement(m_varnames, "y_velocity")) {
345  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 1, mf_comp, 1, 0);
346  mf_comp += 1;
347  }
348  if (containerHasElement(m_varnames, "z_velocity")) {
349  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 2, mf_comp, 1, 0);
350  mf_comp += 1;
351  }
352 
353  if (containerHasElement(m_varnames, "magvel")) {
354 #ifdef _OPENMP
355 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
356 #endif
357  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
358  const amrex::Box& tbx = mfi.tilebox();
359  auto const& dfab = mf_cc_data.array(mfi);
360  auto const& vfab = mf_cc_vel.array(mfi);
361 
362  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
363  {
364  dfab(i,j,k,mf_comp) = std::sqrt(vfab(i,j,k,0)*vfab(i,j,k,0)
365  + vfab(i,j,k,1)*vfab(i,j,k,1)
366  + vfab(i,j,k,2)*vfab(i,j,k,2)) ;
367  });
368  }
369  mf_comp += 1;
370  }
371 
372  if (containerHasElement(m_varnames, "theta")) {
373 #ifdef _OPENMP
374 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
375 #endif
376  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
377  const amrex::Box& tbx = mfi.tilebox();
378  auto const& dfab = mf_cc_data.array(mfi);
379  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
380 
381  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
382  {
383  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoTheta_comp) / cfab(i,j,k,Rho_comp);
384  });
385  }
386  mf_comp += 1;
387  }
388 
389  if (containerHasElement(m_varnames, "qv")) {
390  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp,
391  "qv sampling requested but moisture components not present in state");
392  if (qv_comp >= 0) AMREX_ALWAYS_ASSERT(qv_comp == mf_comp);
393  qv_comp = mf_comp;
394 #ifdef _OPENMP
395 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
396 #endif
397  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
398  const amrex::Box& tbx = mfi.tilebox();
399  auto const& dfab = mf_cc_data.array(mfi);
400  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
401 
402  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
403  {
404  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
405  });
406  }
407  mf_comp += 1;
408  }
409  if (containerHasElement(m_varnames, "qc")) {
410  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ2_comp,
411  "qc sampling requested but moisture components not present in state");
412 #ifdef _OPENMP
413 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
414 #endif
415  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
416  const amrex::Box& tbx = mfi.tilebox();
417  auto const& dfab = mf_cc_data.array(mfi);
418  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
419 
420  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
421  {
422  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ2_comp) / cfab(i,j,k,Rho_comp);
423  });
424  }
425  mf_comp += 1;
426  }
427 
428  if (containerHasElement(m_varnames, "pressure")) {
429  if (vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp) {
430  // moist pressure: use qv from dfab if already sampled, else compute inline
431 #ifdef _OPENMP
432 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
433 #endif
434  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
435  const amrex::Box& tbx = mfi.tilebox();
436  auto const& dfab = mf_cc_data.array(mfi);
437  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
438 
439  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
440  {
441  amrex::Real qv_val = (qv_comp >= 0) ? dfab(i,j,k,qv_comp)
442  : cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
443  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp), qv_val);
444  });
445  }
446  } else {
447  // dry pressure
448 #ifdef _OPENMP
449 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
450 #endif
451  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
452  const amrex::Box& tbx = mfi.tilebox();
453  auto const& dfab = mf_cc_data.array(mfi);
454  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
455 
456  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
457  {
458  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp));
459  });
460  }
461  }
462  mf_comp += 1;
463  }
464 
465  if (containerHasElement(m_varnames, "sgs_tke")) {
466  const bool tke_available =
467  ilev < static_cast<int>(prognostic_tke_available.size()) &&
468  prognostic_tke_available[ilev] != 0;
469  if (!tke_available) {
470  mf_cc_data.setVal(amrex::Real(0.0), mf_comp, 1, 0);
471  } else {
472 #ifdef _OPENMP
473 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
474 #endif
475  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
476  const amrex::Box& tbx = mfi.tilebox();
477  const auto& out = mf_cc_data.array(mfi);
478  const auto& state = vars_new[ilev][Vars::cons].const_array(mfi);
479  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
480  out(i,j,k,mf_comp) = state(i,j,k,RhoKE_comp) / state(i,j,k,Rho_comp);
481  });
482  }
483  }
484  mf_comp += 1;
485  }
486 
487  if (containerHasElement(m_varnames, "sgs_tau13")) {
488  amrex::MultiFab* tau13 = ilev < static_cast<int>(tau13_lev.size())
489  ? tau13_lev[ilev] : nullptr;
490  if (tau13 != nullptr) {
491 #ifdef _OPENMP
492 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
493 #endif
494  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
495  const amrex::Box& tbx = mfi.tilebox();
496  const auto& out = mf_cc_data.array(mfi);
497  const auto& tau = tau13->const_array(mfi);
498  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
499  out(i,j,k,mf_comp) = amrex::Real(0.25) *
500  (tau(i,j,k) + tau(i+1,j,k) +
501  tau(i,j,k+1) + tau(i+1,j,k+1));
502  });
503  }
504  } else {
505  mf_cc_data.setVal(0.0, mf_comp, 1, 0);
506  }
507  mf_comp += 1;
508  }
509 
510  if (containerHasElement(m_varnames, "sgs_tau23")) {
511  amrex::MultiFab* tau23 = ilev < static_cast<int>(tau23_lev.size())
512  ? tau23_lev[ilev] : nullptr;
513  if (tau23 != nullptr) {
514 #ifdef _OPENMP
515 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
516 #endif
517  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
518  const amrex::Box& tbx = mfi.tilebox();
519  const auto& out = mf_cc_data.array(mfi);
520  const auto& tau = tau23->const_array(mfi);
521  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
522  out(i,j,k,mf_comp) = amrex::Real(0.25) *
523  (tau(i,j,k) + tau(i,j+1,k) +
524  tau(i,j,k+1) + tau(i,j+1,k+1));
525  });
526  }
527  } else {
528  mf_cc_data.setVal(0.0, mf_comp, 1, 0);
529  }
530  mf_comp += 1;
531  }
532 
533  if (containerHasElement(m_varnames, "sgs_hfx3")) {
534  amrex::MultiFab* hfx3 = ilev < static_cast<int>(hfx3_lev.size())
535  ? hfx3_lev[ilev] : nullptr;
536  if (hfx3 != nullptr) {
537 #ifdef _OPENMP
538 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
539 #endif
540  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
541  const amrex::Box& tbx = mfi.tilebox();
542  const auto& out = mf_cc_data.array(mfi);
543  const auto& hfx = hfx3->const_array(mfi);
544  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
545  out(i,j,k,mf_comp) = amrex::Real(0.5) * (hfx(i,j,k) + hfx(i,j,k+1));
546  });
547  }
548  } else {
549  mf_cc_data.setVal(0.0, mf_comp, 1, 0);
550  }
551  mf_comp += 1;
552  }
553 
554  m_lev[iline] = ilev;
555  m_ls_mf[iline] = get_line_data(mf_cc_data, dir, cell, bnd_bx);
556 
557  // We can stop if we got the entire line
558  auto min_bnd_bx = m_ls_mf[iline].boxArray().minimalBox();
559  if (bnd_bx == min_bnd_bx) { break; }
560 
561  } // ilev
562 
563  // Fill bnd_bx if use_rbx
564  if (m_use_real_bx) {
565  m_bnd_bx[iline] = m_ls_mf[iline].boxArray().minimalBox();
566  }
567 
568  }// iline
569  }
570 
571  void
572  write_sample_data (amrex::Vector<double>& time,
573  amrex::Vector<int>& level_steps,
574  amrex::Vector<amrex::IntVect>& ref_ratio,
575  amrex::Vector<amrex::Geometry>& geom)
576  {
577  if (m_write_ascii) {
578  write_line_ascii(time);
579  } else {
580  write_line_plotfile(time, level_steps, ref_ratio, geom);
581  }
582  }
583 
584  void
585  write_line_ascii (amrex::Vector<double>& time)
586  {
587  // same as definitions in ERF.H
588  constexpr int datwidth = 14;
589  constexpr int datprecision = 6;
590  constexpr int timeprecision = 13;
591 
592  int nline = static_cast<int>(m_ls_mf.size());
593  int nvar = static_cast<int>(m_varnames.size());
594  for (int iline(0); iline<nline; ++iline) {
595  int dir = m_dir[iline];
596  int lev = m_lev[iline];
597  double m_time = time[lev];
598  amrex::Box m_dom = m_bnd_bx[iline];
599 
600  for (int ivar(0); ivar<nvar; ++ivar) {
601  // Convert multifab to vector
602  amrex::Gpu::HostVector<amrex::Real> vec = sumToLine(m_ls_mf[iline], ivar, 1, m_dom, dir);
603 
604  if (amrex::ParallelDescriptor::IOProcessor()) {
605  int ifile = iline*nvar + ivar;
606  std::ostream& fs = *m_datastream[ifile];
607  fs << std::setw(datwidth) << std::setprecision(timeprecision) << m_time
608  << std::setw(datwidth) << std::setprecision(datprecision);
609  for (const auto& val : vec) {
610  fs << " " << val;
611  }
612  fs << std::endl;
613  }
614  }
615  }
616  }
617 
618  void
619  write_line_plotfile (amrex::Vector<double>& time,
620  amrex::Vector<int>& level_steps,
621  amrex::Vector<amrex::IntVect>& ref_ratio,
622  amrex::Vector<amrex::Geometry>& geom)
623  {
624  int nline = static_cast<int>(m_ls_mf.size());
625  for (int iline(0); iline<nline; ++iline) {
626  // Data members that can be used as-is
627  int dir = m_dir[iline];
628  int lev = m_lev[iline];
629  double m_time = time[lev];
630  amrex::Vector<int> m_level_steps = {level_steps[lev]};
631  amrex::Vector<amrex::IntVect> m_ref_ratio = {ref_ratio[lev]};
632 
633  // Create modified geometry object corresponding to the line
634  auto plo = geom[lev].ProbLo();
635  auto dx = geom[lev].CellSize();
636  amrex::Vector<amrex::Geometry> m_geom; m_geom.resize(1);
637  amrex::Vector<int> is_per(AMREX_SPACEDIM,0);
638  amrex::Box m_dom = m_bnd_bx[iline];
639  amrex::RealBox m_rb;
640  for (int d(0); d<AMREX_SPACEDIM; ++d) {
641  amrex::Real offset = (d==dir) ? 0 : myhalf;
642  amrex::Real lo = plo[d] + ( m_dom.smallEnd(d) - offset ) * dx[d];
643  amrex::Real hi = plo[d] + ( m_dom.bigEnd(d) + offset ) * dx[d];
644 
645  m_rb.setLo(d,lo);
646  m_rb.setHi(d,hi);
647 
648  is_per[d] = geom[lev].isPeriodic(d);
649  }
650  m_geom[0].define(m_dom, &m_rb, geom[lev].Coord(), is_per.data());
651 
652  // Create plotfile name
653  std::string name_line = m_name[iline];
654  name_line += "_step_";
655  std::string plotfilename = amrex::Concatenate(name_line, m_level_steps[0], 5);
656 
657  // Get the data
658  amrex::Vector<const amrex::MultiFab*> mf = {&(m_ls_mf[iline])};
659 
660  // Write each line
661  WriteMultiLevelPlotfile(plotfilename, 1, mf,
662  m_varnames, m_geom, static_cast<amrex::Real>(m_time),
663  m_level_steps, m_ref_ratio);
664  }
665  }
666 
667  amrex::Vector<int> m_dir;
668  amrex::Vector<int> m_lev;
669  amrex::Vector<amrex::Box> m_bnd_bx;
670  amrex::Vector<amrex::RealBox> m_bnd_rbx;
671  amrex::Vector<amrex::MultiFab> m_ls_mf;
672  amrex::Vector<std::string> m_name;
673 
674  bool m_use_real_bx{false};
675  bool m_write_ascii{false};
676  amrex::Vector<std::string> m_varnames {"magvel","theta"};
677  amrex::Vector<std::unique_ptr<std::fstream> > m_datastream;
678 };
679 
680 
681 struct PlaneSampler
682 {
683  PlaneSampler ()
684  {
685  amrex::ParmParse pp("erf");
686 
687  // Count number of lo and hi points define the plane
688  int n_plane_lo = pp.countval("sample_plane_lo") / AMREX_SPACEDIM;
689  int n_plane_hi = pp.countval("sample_plane_hi") / AMREX_SPACEDIM;
690  int n_plane_dir = pp.countval("sample_plane_dir");
691  AMREX_ALWAYS_ASSERT( (n_plane_lo==n_plane_hi ) &&
692  (n_plane_lo==n_plane_dir) );
693 
694  // Parse the data
695  if (n_plane_lo > 0) {
696  // Parse lo
697  amrex::Vector<amrex::Real> r_lo; r_lo.resize(n_plane_lo*AMREX_SPACEDIM);
698  amrex::Vector<amrex::Vector<amrex::Real>> rv_lo;
699  pp.queryarr("sample_plane_lo",r_lo,0,n_plane_lo*AMREX_SPACEDIM);
700  for (int i(0); i < n_plane_lo; i++) {
701  amrex::Vector<amrex::Real> rv = {r_lo[AMREX_SPACEDIM*i+0],
702  r_lo[AMREX_SPACEDIM*i+1],
703  r_lo[AMREX_SPACEDIM*i+2]};
704  rv_lo.push_back(rv);
705  }
706 
707  // Parse hi
708  amrex::Vector<amrex::Real> r_hi; r_hi.resize(n_plane_hi*AMREX_SPACEDIM);
709  amrex::Vector<amrex::Vector<amrex::Real>> rv_hi;
710  pp.queryarr("sample_plane_hi",r_hi,0,n_plane_hi*AMREX_SPACEDIM);
711  for (int i(0); i < n_plane_hi; i++) {
712  amrex::Vector<amrex::Real> rv = {r_hi[AMREX_SPACEDIM*i+0],
713  r_hi[AMREX_SPACEDIM*i+1],
714  r_hi[AMREX_SPACEDIM*i+2]};
715  rv_hi.push_back(rv);
716  }
717 
718  // Construct vector of bounding real boxes
719  m_bnd_rbx.resize(n_plane_lo);
720  for (int i(0); i < n_plane_hi; i++){
721  amrex::RealBox rbx(rv_lo[i].data(),rv_hi[i].data());
722  m_bnd_rbx[i] = rbx;
723  }
724 
725  // Parse directionality
726  m_dir.resize(n_plane_dir);
727  pp.queryarr("sample_plane_dir",m_dir,0,n_plane_dir);
728 
729  // Parse names
730  std::string name_base = "plt_plane_";
731  m_name.resize(n_plane_lo);
732  int n_names = pp.countval("sample_plane_name");
733  if (n_names > 0) {
734  AMREX_ALWAYS_ASSERT( n_names==n_plane_lo );
735  pp.queryarr("sample_plane_name",m_name,0,n_names);
736  } else {
737  for (int iplane(0); iplane<n_plane_lo; ++iplane) {
738  m_name[iplane] = amrex::Concatenate(name_base, iplane , 5);
739  }
740  }
741 
742  // Optional cap on the finest level written (-1 => all intersecting levels)
743  pp.queryAdd("plane_sampling_max_level", m_max_level);
744 
745  // Allocate per-plane vectors of per-level MF pointers
746  m_ps_mf.resize(n_plane_lo);
747 
748  // Get requested vars
749  if (pp.countval("plane_sampling_vars") > 0) {
750  m_varnames.clear();
751  amrex::Vector<std::string> requested_vars;
752  pp.queryarr("plane_sampling_vars",requested_vars);
753  amrex::Print() << "Selected plane sampling vars :";
754  if (containerHasElement(requested_vars, "density")) {
755  m_varnames.push_back("density");
756  amrex::Print() << " " << "density";
757  }
758  if (containerHasElement(requested_vars, "x_velocity")) {
759  m_varnames.push_back("x_velocity");
760  amrex::Print() << " " << "x_velocity";
761  }
762  if (containerHasElement(requested_vars, "y_velocity")) {
763  m_varnames.push_back("y_velocity");
764  amrex::Print() << " " << "y_velocity";
765  }
766  if (containerHasElement(requested_vars, "z_velocity")) {
767  m_varnames.push_back("z_velocity");
768  amrex::Print() << " " << "z_velocity";
769  }
770  if (containerHasElement(requested_vars, "magvel")) {
771  m_varnames.push_back("magvel");
772  amrex::Print() << " " << "magvel";
773  }
774  if (containerHasElement(requested_vars, "theta")) {
775  m_varnames.push_back("theta");
776  amrex::Print() << " " << "theta";
777  }
778  if (containerHasElement(requested_vars, "qv")) {
779  m_varnames.push_back("qv");
780  amrex::Print() << " " << "qv";
781  }
782  if (containerHasElement(requested_vars, "qc")) {
783  m_varnames.push_back("qc");
784  amrex::Print() << " " << "qc";
785  }
786  if (containerHasElement(requested_vars, "pressure")) {
787  m_varnames.push_back("pressure");
788  amrex::Print() << " " << "pressure";
789  }
790  amrex::Print() << std::endl;
791  }
792  }
793  }
794 
795  // This must match what is in AMReX_MultiFabUtil.H
796  amrex::Box
797  getIndexBox (const amrex::RealBox& real_box,
798  const amrex::Geometry& geom) {
799  amrex::IntVect slice_lo, slice_hi;
800 
801  AMREX_D_TERM(slice_lo[0]=static_cast<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
802  slice_lo[1]=static_cast<int>(std::floor((real_box.lo(1) - geom.ProbLo(1))/geom.CellSize(1)));,
803  slice_lo[2]=static_cast<int>(std::floor((real_box.lo(2) - geom.ProbLo(2))/geom.CellSize(2))););
804 
805  AMREX_D_TERM(slice_hi[0]=static_cast<int>(std::floor((real_box.hi(0) - geom.ProbLo(0))/geom.CellSize(0)));,
806  slice_hi[1]=static_cast<int>(std::floor((real_box.hi(1) - geom.ProbLo(1))/geom.CellSize(1)));,
807  slice_hi[2]=static_cast<int>(std::floor((real_box.hi(2) - geom.ProbLo(2))/geom.CellSize(2))););
808 
809  return amrex::Box(slice_lo, slice_hi) & geom.Domain();
810  }
811 
812  void
813  get_sample_data (amrex::Vector<amrex::Geometry>& geom,
814  amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new)
815  {
816  int nlev = static_cast<int>(vars_new.size());
817  int nplane = static_cast<int>(m_bnd_rbx.size());
818  int ncomp = static_cast<int>(m_varnames.size());
819  bool interpolate = true;
820 
821  int qv_comp = -1;
822 
823  // Loop over each plane
824  for (int iplane(0); iplane<nplane; ++iplane) {
825  int dir = m_dir[iplane];
826  amrex::RealBox bnd_rbx = m_bnd_rbx[iplane];
827  amrex::Real point = bnd_rbx.lo(dir);
828 
829  // Collect every level 0..lev_cap whose grids intersect the plane.
830  int lev_cap = (m_max_level >= 0 && m_max_level < nlev-1) ? m_max_level : nlev-1;
831  m_ps_mf[iplane].clear();
832 
833  for (int ilev(0); ilev<=lev_cap; ++ilev) {
834 
835  // Stop ascending once a level's grids no longer touch the flattened plane
836  {
837  amrex::Box plane_bx = getIndexBox(bnd_rbx, geom[ilev]);
838  int k_l = static_cast<int>(std::floor((point - geom[ilev].ProbLo(dir))
839  / geom[ilev].CellSize(dir)));
840  plane_bx.setSmall(dir, k_l); plane_bx.setBig(dir, k_l);
841  if (!vars_new[ilev][Vars::cons].boxArray().intersects(plane_bx)) { break; }
842  }
843 
844  // Construct CC velocities
845  amrex::MultiFab mf_cc_vel;
846  auto ba = vars_new[ilev][Vars::cons].boxArray();
847  auto dm = vars_new[ilev][Vars::cons].DistributionMap();
848  mf_cc_vel.define(ba, dm, AMREX_SPACEDIM, amrex::IntVect(1,1,1));
849  average_face_to_cellcenter(mf_cc_vel,0,
850  amrex::Array<const amrex::MultiFab*,3>{&vars_new[ilev][Vars::xvel],
851  &vars_new[ilev][Vars::yvel],
852  &vars_new[ilev][Vars::zvel]});
853 
854  // Construct MultiFab holding requested variables
855  amrex::MultiFab mf_cc_data;
856  mf_cc_data.define(ba, dm, ncomp, 1);
857 
858  int mf_comp = 0;
859 
860  if (containerHasElement(m_varnames, "density")) {
861  amrex::MultiFab::Copy(mf_cc_data, vars_new[ilev][Vars::cons], Rho_comp, mf_comp, 1, 0);
862  mf_comp += 1;
863  }
864 
865  if (containerHasElement(m_varnames, "x_velocity")) {
866  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 0, mf_comp, 1, 0);
867  mf_comp += 1;
868  }
869  if (containerHasElement(m_varnames, "y_velocity")) {
870  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 1, mf_comp, 1, 0);
871  mf_comp += 1;
872  }
873  if (containerHasElement(m_varnames, "z_velocity")) {
874  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 2, mf_comp, 1, 0);
875  mf_comp += 1;
876  }
877 
878  if (containerHasElement(m_varnames, "magvel")) {
879 #ifdef _OPENMP
880 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
881 #endif
882  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
883  const amrex::Box& tbx = mfi.tilebox();
884  auto const& dfab = mf_cc_data.array(mfi);
885  auto const& vfab = mf_cc_vel.array(mfi);
886 
887  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
888  {
889  dfab(i,j,k,mf_comp) = std::sqrt(vfab(i,j,k,0)*vfab(i,j,k,0)
890  + vfab(i,j,k,1)*vfab(i,j,k,1)
891  + vfab(i,j,k,2)*vfab(i,j,k,2));
892  });
893  }
894  mf_comp += 1;
895  }
896 
897  if (containerHasElement(m_varnames, "theta")) {
898 #ifdef _OPENMP
899 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
900 #endif
901  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
902  const amrex::Box& tbx = mfi.tilebox();
903  auto const& dfab = mf_cc_data.array(mfi);
904  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
905 
906  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
907  {
908  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoTheta_comp) / cfab(i,j,k,Rho_comp);
909  });
910  }
911  mf_comp += 1;
912  }
913 
914  if (containerHasElement(m_varnames, "qv")) {
915  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp,
916  "qv sampling requested but moisture components not present in state");
917  if (qv_comp >= 0) AMREX_ALWAYS_ASSERT(qv_comp == mf_comp);
918  qv_comp = mf_comp;
919 #ifdef _OPENMP
920 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
921 #endif
922  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
923  const amrex::Box& tbx = mfi.tilebox();
924  auto const& dfab = mf_cc_data.array(mfi);
925  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
926 
927  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
928  {
929  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
930  });
931  }
932  mf_comp += 1;
933  }
934  if (containerHasElement(m_varnames, "qc")) {
935  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ2_comp,
936  "qc sampling requested but moisture components not present in state");
937 #ifdef _OPENMP
938 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
939 #endif
940  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
941  const amrex::Box& tbx = mfi.tilebox();
942  auto const& dfab = mf_cc_data.array(mfi);
943  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
944 
945  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
946  {
947  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ2_comp) / cfab(i,j,k,Rho_comp);
948  });
949  }
950  mf_comp += 1;
951  }
952 
953  if (containerHasElement(m_varnames, "pressure")) {
954  if (vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp) {
955  // moist pressure: use qv from dfab if already sampled, else compute inline
956 #ifdef _OPENMP
957 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
958 #endif
959  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
960  const amrex::Box& tbx = mfi.tilebox();
961  auto const& dfab = mf_cc_data.array(mfi);
962  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
963 
964  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
965  {
966  amrex::Real qv_val = (qv_comp >= 0) ? dfab(i,j,k,qv_comp)
967  : cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
968  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp), qv_val);
969  });
970  }
971  } else {
972  // dry pressure
973 #ifdef _OPENMP
974 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
975 #endif
976  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
977  const amrex::Box& tbx = mfi.tilebox();
978  auto const& dfab = mf_cc_data.array(mfi);
979  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
980 
981  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
982  {
983  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp));
984  });
985  }
986  }
987  mf_comp += 1;
988  }
989 
990  auto slice = get_slice_data(dir, point, mf_cc_data, geom[ilev],
991  0, ncomp, interpolate, bnd_rbx);
992 
993  // Defensive: an empty return means nothing was actually sampled.
994  if (!slice || slice->boxArray().size() == 0) { break; }
995 
996  // Broadcast the single sampled plane across 2^ilev cells in dir so
997  // every level shares an isotropic refinement ratio. The dir extent
998  // becomes [0, 2^ilev-1]; each cell holds the interpolated plane value
999  // replicated across the shared level-0 slab thickness.
1000  const int rr_l = 1 << ilev;
1001  const int k_slice = slice->boxArray().minimalBox().smallEnd(dir);
1002  const int l_dir = dir;
1003 
1004  amrex::BoxList bl;
1005  const amrex::BoxArray& slice_ba = slice->boxArray();
1006  for (int ib(0); ib<static_cast<int>(slice_ba.size()); ++ib) {
1007  amrex::Box b = slice_ba[ib];
1008  b.setSmall(dir, 0);
1009  b.setBig(dir, rr_l - 1);
1010  bl.push_back(b);
1011  }
1012  amrex::BoxArray out_ba(std::move(bl));
1013  auto out_mf = std::make_unique<amrex::MultiFab>(out_ba, slice->DistributionMap(),
1014  ncomp, 0);
1015 
1016  for (amrex::MFIter mfi(*out_mf); mfi.isValid(); ++mfi) {
1017  const amrex::Box& obx = mfi.validbox();
1018  auto const& ofab = out_mf->array(mfi);
1019  auto const& sfab = slice->array(mfi);
1020  amrex::ParallelFor(obx, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept
1021  {
1022  int si = i, sj = j, sk = k;
1023  if (l_dir == 0) { si = k_slice; }
1024  else if (l_dir == 1) { sj = k_slice; }
1025  else { sk = k_slice; }
1026  ofab(i,j,k,n) = sfab(si,sj,sk,n);
1027  });
1028  }
1029 
1030  m_ps_mf[iplane].push_back(std::move(out_mf));
1031 
1032  } // ilev
1033  }// iplane
1034  }
1035 
1036  void
1037  write_sample_data (amrex::Vector<double>& time,
1038  amrex::Vector<int>& level_steps,
1039  amrex::Vector<amrex::IntVect>& ref_ratio,
1040  amrex::Vector<amrex::Geometry>& geom)
1041  {
1042  amrex::ignore_unused(ref_ratio);
1043  int nplane = m_ps_mf.size();
1044  for (int iplane(0); iplane<nplane; ++iplane) {
1045  int nlev_c = static_cast<int>(m_ps_mf[iplane].size());
1046  if (nlev_c == 0) { continue; }
1047 
1048  int dir = m_dir[iplane];
1049  amrex::RealBox bnd_rbx = m_bnd_rbx[iplane];
1050  amrex::Real point = bnd_rbx.lo(dir);
1051 
1052  // One shared physical RealBox for every level: full in-plane extent of the
1053  // requested plane, and a dir-slab of level-0 thickness centered on the plane.
1054  // Level 0 supplies prob_lo/hi to the writer, so this object must be shared.
1055  amrex::Real dx0 = geom[0].CellSize(dir);
1056  amrex::RealBox shared_rbx = bnd_rbx;
1057  shared_rbx.setLo(dir, point - myhalf*dx0);
1058  shared_rbx.setHi(dir, point + myhalf*dx0);
1059 
1060  amrex::Vector<int> is_per(AMREX_SPACEDIM,0);
1061  for (int d(0); d<AMREX_SPACEDIM; ++d) { is_per[d] = geom[0].isPeriodic(d); }
1062 
1063  // Build the index chain by refining the level-0 plane box so that
1064  // Domain_l == refine(Domain_{l-1},2) holds exactly in every axis.
1065  amrex::Box B0 = getIndexBox(bnd_rbx, geom[0]);
1066 
1067  amrex::Vector<const amrex::MultiFab*> mf(nlev_c);
1068  amrex::Vector<amrex::Geometry> m_geom(nlev_c);
1069  amrex::Vector<int> m_level_steps(nlev_c);
1070  amrex::Vector<amrex::IntVect> m_ref_ratio(nlev_c-1, amrex::IntVect(2));
1071 
1072  for (int l(0); l<nlev_c; ++l) {
1073  const int rr_l = 1 << l;
1074 
1075  // The nested plane plotfile stores a single scalar refinement
1076  // ratio per level, so multi-level output requires factor-2
1077  // isotropic refinement.
1078  for (int d(0); d<AMREX_SPACEDIM; ++d) {
1079  amrex::Real ratio = geom[0].CellSize(d) / geom[l].CellSize(d);
1080  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
1081  std::abs(ratio - static_cast<amrex::Real>(rr_l)) < amrex::Real(1.e-6) * rr_l,
1082  "PlaneSampler multi-level output requires factor-2 isotropic refinement (amr.ref_ratio = 2 2 2)");
1083  }
1084 
1085  amrex::Box domain_l = amrex::refine(B0, rr_l);
1086  domain_l.setSmall(dir, 0);
1087  domain_l.setBig(dir, rr_l - 1);
1088 
1089  mf[l] = m_ps_mf[iplane][l].get();
1090  m_geom[l].define(domain_l, &shared_rbx, geom[l].Coord(), is_per.data());
1091  m_level_steps[l] = level_steps[l];
1092  AMREX_ASSERT(domain_l.contains(mf[l]->boxArray().minimalBox()));
1093  }
1094 
1095  // Create plotfile name
1096  std::string name_plane = m_name[iplane];
1097  name_plane += "_step_";
1098  std::string plotfilename = amrex::Concatenate(name_plane, m_level_steps[0], 5);
1099 
1100  // Write the nested multi-level plane plotfile
1101  WriteMultiLevelPlotfile(plotfilename, nlev_c, mf,
1102  m_varnames, m_geom, static_cast<amrex::Real>(time[0]),
1103  m_level_steps, m_ref_ratio);
1104  } // iplane
1105  }
1106 
1107  int m_max_level = -1;
1108  amrex::Vector<int> m_dir;
1109  amrex::Vector<amrex::RealBox> m_bnd_rbx;
1110  amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> m_ps_mf;
1111  amrex::Vector<std::string> m_name;
1112 
1113  amrex::Vector<std::string> m_varnames {"magvel","theta"};
1114 };
1115 #endif
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
ERF_Container.H
containerHasElement
bool containerHasElement(const V &iterable, const T &query)
Definition: ERF_Container.H:5
tau23
@ tau23
Definition: ERF_DataStruct.H:39
tau13
@ tau13
Definition: ERF_DataStruct.H:39
Coord
Coord
Coordinate-axis selector.
Definition: ERF_DataStruct.H:144
ERF_EOS.H
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
ERF_IndexDefines.H
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:46
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:45
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:41
pp
ParmParse pp("prob")
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_cloud_chamber_config.active, "Cloud Chamber: initializer reached without a parsed configuration")
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:31
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
MORRInd::fs
@ fs
Definition: ERF_AdvanceMorrison.cpp:122
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:217
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216
LineSampler
Definition: ERF_SampleData.H:21
LineSampler::write_line_plotfile
void write_line_plotfile(amrex::Vector< double > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:619
LineSampler::m_varnames
amrex::Vector< std::string > m_varnames
Definition: ERF_SampleData.H:676
LineSampler::m_ls_mf
amrex::Vector< amrex::MultiFab > m_ls_mf
Definition: ERF_SampleData.H:671
LineSampler::get_sample_data
void get_sample_data(amrex::Vector< amrex::Geometry > &geom, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_new, const amrex::Vector< amrex::MultiFab * > &tau13_lev={}, const amrex::Vector< amrex::MultiFab * > &tau23_lev={}, const amrex::Vector< amrex::MultiFab * > &hfx3_lev={}, const amrex::Vector< int > &prognostic_tke_available={})
Definition: ERF_SampleData.H:295
LineSampler::write_coords
void write_coords(amrex::Vector< std::unique_ptr< amrex::MultiFab > > &z_phys_cc, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:242
LineSampler::write_sample_data
void write_sample_data(amrex::Vector< double > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:572
LineSampler::m_datastream
amrex::Vector< std::unique_ptr< std::fstream > > m_datastream
Definition: ERF_SampleData.H:677
LineSampler::m_use_real_bx
bool m_use_real_bx
Definition: ERF_SampleData.H:674
LineSampler::write_line_ascii
void write_line_ascii(amrex::Vector< double > &time)
Definition: ERF_SampleData.H:585
LineSampler::m_dir
amrex::Vector< int > m_dir
Definition: ERF_SampleData.H:667
LineSampler::m_name
amrex::Vector< std::string > m_name
Definition: ERF_SampleData.H:672
LineSampler::LineSampler
LineSampler()
Definition: ERF_SampleData.H:22
LineSampler::m_bnd_bx
amrex::Vector< amrex::Box > m_bnd_bx
Definition: ERF_SampleData.H:669
LineSampler::getIndexBox
amrex::Box getIndexBox(const amrex::RealBox &real_box, const amrex::Geometry &geom)
Definition: ERF_SampleData.H:226
LineSampler::m_write_ascii
bool m_write_ascii
Definition: ERF_SampleData.H:675
LineSampler::m_bnd_rbx
amrex::Vector< amrex::RealBox > m_bnd_rbx
Definition: ERF_SampleData.H:670
LineSampler::m_lev
amrex::Vector< int > m_lev
Definition: ERF_SampleData.H:668
PlaneSampler
Definition: ERF_SampleData.H:682
PlaneSampler::m_dir
amrex::Vector< int > m_dir
Definition: ERF_SampleData.H:1108
PlaneSampler::getIndexBox
amrex::Box getIndexBox(const amrex::RealBox &real_box, const amrex::Geometry &geom)
Definition: ERF_SampleData.H:797
PlaneSampler::write_sample_data
void write_sample_data(amrex::Vector< double > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:1037
PlaneSampler::m_max_level
int m_max_level
Definition: ERF_SampleData.H:1107
PlaneSampler::m_name
amrex::Vector< std::string > m_name
Definition: ERF_SampleData.H:1111
PlaneSampler::m_ps_mf
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_ps_mf
Definition: ERF_SampleData.H:1110
PlaneSampler::m_bnd_rbx
amrex::Vector< amrex::RealBox > m_bnd_rbx
Definition: ERF_SampleData.H:1109
PlaneSampler::PlaneSampler
PlaneSampler()
Definition: ERF_SampleData.H:683
PlaneSampler::m_varnames
amrex::Vector< std::string > m_varnames
Definition: ERF_SampleData.H:1113
PlaneSampler::get_sample_data
void get_sample_data(amrex::Vector< amrex::Geometry > &geom, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_new)
Definition: ERF_SampleData.H:813