ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_SampleData.H
Go to the documentation of this file.
1 #ifndef ERF_SAMPLEDATA_H
2 #define ERF_SAMPLEDATA_H
3 
4 #include <memory>
5 
6 #include <AMReX_ParmParse.H>
7 #include <AMReX_MultiFab.H>
8 #include <AMReX_MultiFabUtil.H>
9 #include <AMReX_PlotFileUtil.H>
10 
11 #include <ERF_IndexDefines.H>
12 #include <ERF_EOS.H>
13 #include <ERF_Container.H>
14 
15 #include <iostream>
16 
18 {
20  {
21  amrex::ParmParse pp("erf");
22 
23  bool has_line_idx = pp.contains("sample_line_lo") || pp.contains("sample_line_hi");
24  bool has_line_real = pp.contains("sample_line_lo_real") || pp.contains("sample_line_hi_real");
25  if (has_line_idx && has_line_real) {
26  amrex::Abort("Specify only one of erf.sample_line_lo/hi or erf.sample_line_lo_real/hi_real");
27  }
28 
29  // Count number of lo and hi points define the line
30  int n_line_lo = pp.countval("sample_line_lo") / AMREX_SPACEDIM;
31  int n_line_hi = pp.countval("sample_line_hi") / AMREX_SPACEDIM;
32  int n_line_lo_real = pp.countval("sample_line_lo_real") / AMREX_SPACEDIM;
33  int n_line_hi_real = pp.countval("sample_line_hi_real") / AMREX_SPACEDIM;
34  int n_line_dir = pp.countval("sample_line_dir");
35  if (has_line_idx) {
36  AMREX_ALWAYS_ASSERT( (n_line_lo==n_line_hi ) &&
37  (n_line_lo==n_line_dir) );
38  } else if (has_line_real) {
39  AMREX_ALWAYS_ASSERT( (n_line_lo_real==n_line_hi_real ) &&
40  (n_line_lo_real==n_line_dir) );
41  }
42 
43  // Parse the data
44  int nline = 0;
45  if (has_line_idx) {
46  nline = n_line_lo;
47  if (nline > 0) {
48  // Parse lo
49  amrex::Vector<int> idx_lo; idx_lo.resize(nline*AMREX_SPACEDIM);
50  amrex::Vector<amrex::IntVect> iv_lo; iv_lo.resize(nline);
51  pp.queryarr("sample_line_lo",idx_lo,0,nline*AMREX_SPACEDIM);
52  for (int i(0); i < nline; i++) {
53  amrex::IntVect iv(idx_lo[AMREX_SPACEDIM*i+0],
54  idx_lo[AMREX_SPACEDIM*i+1],
55  idx_lo[AMREX_SPACEDIM*i+2]);
56  iv_lo[i] = iv;
57  }
58 
59  // Parse hi
60  amrex::Vector<int> idx_hi; idx_hi.resize(nline*AMREX_SPACEDIM);
61  amrex::Vector<amrex::IntVect> iv_hi; iv_hi.resize(nline);
62  pp.queryarr("sample_line_hi",idx_hi,0,nline*AMREX_SPACEDIM);
63  for (int i(0); i < nline; i++) {
64  amrex::IntVect iv(idx_hi[AMREX_SPACEDIM*i+0],
65  idx_hi[AMREX_SPACEDIM*i+1],
66  idx_hi[AMREX_SPACEDIM*i+2]);
67  iv_hi[i] = iv;
68  }
69 
70  // Construct vector of bounding boxes
71  m_bnd_bx.resize(nline);
72  for (int i = 0; i < nline; i++){
73  amrex::Box lbx(iv_lo[i],iv_hi[i]);
74  m_bnd_bx[i] = lbx;
75  }
76  }
77  } else if (has_line_real) {
78  m_use_real_bx = true;
79  nline = n_line_lo_real;
80  if (nline > 0) {
81  // Parse lo
82  amrex::Vector<amrex::Real> real_lo; real_lo.resize(nline*AMREX_SPACEDIM);
83  amrex::Vector<amrex::Vector<amrex::Real>> rv_lo;
84  pp.queryarr("sample_line_lo_real",real_lo,0,nline*AMREX_SPACEDIM);
85  for (int i(0); i < nline; i++) {
86  amrex::Vector<amrex::Real> rv = {real_lo[AMREX_SPACEDIM*i+0],
87  real_lo[AMREX_SPACEDIM*i+1],
88  real_lo[AMREX_SPACEDIM*i+2]};
89  rv_lo.push_back(rv);
90  }
91 
92  // Parse hi
93  amrex::Vector<amrex::Real> real_hi; real_hi.resize(nline*AMREX_SPACEDIM);
94  amrex::Vector<amrex::Vector<amrex::Real>> rv_hi;
95  pp.queryarr("sample_line_hi_real",real_hi,0,nline*AMREX_SPACEDIM);
96  for (int i(0); i < nline; i++) {
97  amrex::Vector<amrex::Real> rv = {real_hi[AMREX_SPACEDIM*i+0],
98  real_hi[AMREX_SPACEDIM*i+1],
99  real_hi[AMREX_SPACEDIM*i+2]};
100  rv_hi.push_back(rv);
101  }
102 
103  // Construct vector of bounding real boxes
104  m_bnd_rbx.resize(nline);
105  for (int i = 0; i < nline; i++){
106  amrex::RealBox rbx(rv_lo[i].data(),rv_hi[i].data());
107  m_bnd_rbx[i] = rbx;
108  }
109 
110  // Construct vector of bounding boxes
111  m_bnd_bx.resize(nline);
112  }
113  }
114 
115  if (nline > 0) {
116  // Parse directionality
117  m_dir.resize(n_line_dir);
118  pp.queryarr("sample_line_dir",m_dir,0,n_line_dir);
119 
120  // Parse names
121  std::string name_base = "plt_line_";
122  m_name.resize(nline);
123  int n_names = pp.countval("sample_line_name");
124  if (n_names > 0) {
125  AMREX_ALWAYS_ASSERT( n_names==nline );
126  pp.queryarr("sample_line_name",m_name,0,n_names);
127  } else {
128  for (int iline(0); iline<nline; ++iline) {
129  m_name[iline] = amrex::Concatenate(name_base, iline , 5);
130  }
131  }
132 
133  // Allocate space for level indicator
134  m_lev.resize(n_line_dir,0);
135 
136  // Allocate space for MF pointers
137  m_ls_mf.resize(nline);
138 
139  // Get requested vars
140  if (pp.countval("line_sampling_vars") > 0) {
141  m_varnames.clear();
142  amrex::Vector<std::string> requested_vars;
143  pp.queryarr("line_sampling_vars",requested_vars);
144  amrex::Print() << "Selected line sampling vars :";
145  if (containerHasElement(requested_vars, "density")) {
146  m_varnames.push_back("density");
147  amrex::Print() << " " << "density";
148  }
149  if (containerHasElement(requested_vars, "x_velocity")) {
150  m_varnames.push_back("x_velocity");
151  amrex::Print() << " " << "x_velocity";
152  }
153  if (containerHasElement(requested_vars, "y_velocity")) {
154  m_varnames.push_back("y_velocity");
155  amrex::Print() << " " << "y_velocity";
156  }
157  if (containerHasElement(requested_vars, "z_velocity")) {
158  m_varnames.push_back("z_velocity");
159  amrex::Print() << " " << "z_velocity";
160  }
161  if (containerHasElement(requested_vars, "magvel")) {
162  m_varnames.push_back("magvel");
163  amrex::Print() << " " << "magvel";
164  }
165  if (containerHasElement(requested_vars, "theta")) {
166  m_varnames.push_back("theta");
167  amrex::Print() << " " << "theta";
168  }
169  if (containerHasElement(requested_vars, "qv")) {
170  m_varnames.push_back("qv");
171  amrex::Print() << " " << "qv";
172  }
173  if (containerHasElement(requested_vars, "qc")) {
174  m_varnames.push_back("qc");
175  amrex::Print() << " " << "qc";
176  }
177  if (containerHasElement(requested_vars, "pressure")) {
178  m_varnames.push_back("pressure");
179  amrex::Print() << " " << "pressure";
180  }
181  amrex::Print() << std::endl;
182  }
183 
184  // Write outputs to text files, one file per variable, with all
185  // times appended to the same file
186  pp.query("line_sampling_text_output",m_write_ascii);
187  if (m_write_ascii && amrex::ParallelDescriptor::IOProcessor()) {
188  int nvar = static_cast<int>(m_varnames.size());
189  m_datastream.resize(nline * nvar);
190  int i = 0;
191  for (int iline(0); iline<nline; ++iline) {
192  for (int ivar(0); ivar<nvar; ++ivar) {
193  std::string filename = m_name[iline] + "." + m_varnames[ivar];
194  m_datastream[i] = std::make_unique<std::fstream>();
195  m_datastream[i]->open(filename.c_str(),std::ios::out|std::ios::app);
196  if (!m_datastream[i]->good()) {
197  amrex::FileOpenFailed(filename);
198  }
199  i++;
200  }
201  }
202  }
203  }
204  }
205 
206  amrex::Box
207  getIndexBox (const amrex::RealBox& real_box,
208  const amrex::Geometry& geom) {
209  amrex::IntVect slice_lo, slice_hi;
210 
211  AMREX_D_TERM(slice_lo[0]=static_cast<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
212  slice_lo[1]=static_cast<int>(std::floor((real_box.lo(1) - geom.ProbLo(1))/geom.CellSize(1)));,
213  slice_lo[2]=static_cast<int>(std::floor((real_box.lo(2) - geom.ProbLo(2))/geom.CellSize(2))););
214 
215  AMREX_D_TERM(slice_hi[0]=static_cast<int>(std::floor((real_box.hi(0) - geom.ProbLo(0))/geom.CellSize(0)));,
216  slice_hi[1]=static_cast<int>(std::floor((real_box.hi(1) - geom.ProbLo(1))/geom.CellSize(1)));,
217  slice_hi[2]=static_cast<int>(std::floor((real_box.hi(2) - geom.ProbLo(2))/geom.CellSize(2))););
218 
219  return amrex::Box(slice_lo, slice_hi) & geom.Domain();
220  }
221 
222  void
223  write_coords (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& z_phys_cc,
224  amrex::Vector<amrex::Geometry>& geom)
225  {
226  if (!m_write_ascii) return;
227 
228  amrex::Print() << "Writing out line coordinates to text" << std::endl;
229 
230  for (int lev(0); lev < z_phys_cc.size(); ++lev) {
231  // Write one text file per level
232  std::ofstream outfile;
233  if (amrex::ParallelDescriptor::IOProcessor()) {
234  std::string fname = amrex::Concatenate("plt_line_lev", lev, 1);
235  fname += ".zcc";
236  outfile.open(fname);
237 
238  if (!outfile.is_open()) {
239  amrex::AllPrint() << "Could not open " << fname << std::endl;
240  }
241  }
242 
243  // Loop over each line
244  int nline = static_cast<int>(m_ls_mf.size());
245  for (int iline(0); iline<nline; ++iline) {
246  int dir = m_dir[iline];
247  amrex::Box bnd_bx = (m_use_real_bx)
248  ? getIndexBox(m_bnd_rbx[iline], geom[lev])
249  : m_bnd_bx[iline];
250  amrex::IntVect first_cell = bnd_bx.smallEnd();
251 
252  // Create multifab with "sampled" z_phys values
253  amrex::MultiFab line_coords_mf = get_line_data(
254  *z_phys_cc[lev], dir, first_cell, bnd_bx
255  );
256 
257  // Convert multifab to vector
258  amrex::Gpu::HostVector<amrex::Real> vec = sumToLine(
259  line_coords_mf, 0, 1, bnd_bx, dir
260  );
261 
262  // Append to file
263  if (amrex::ParallelDescriptor::IOProcessor()) {
264  for (const auto& zval : vec) {
265  outfile << " " << zval;
266  }
267  outfile << std::endl;
268  }
269  } // line loop
270 
271  outfile.close();
272  } // level loop
273  }
274 
275  void
276  get_sample_data (amrex::Vector<amrex::Geometry>& geom,
277  amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new)
278  {
279  int nlev = static_cast<int>(vars_new.size());
280  int nline = static_cast<int>(m_bnd_bx.size());
281  int ncomp = static_cast<int>(m_varnames.size());
282 
283  int qv_comp = -1;
284 
285  // Loop over each line
286  for (int iline(0); iline<nline; ++iline) {
287  int dir = m_dir[iline];
288 
289  // Search each level to get the finest data possible
290  for (int ilev(nlev-1); ilev>=0; --ilev) {
291  amrex::Box bnd_bx = (m_use_real_bx)
292  ? getIndexBox(m_bnd_rbx[iline], geom[ilev])
293  : m_bnd_bx[iline];
294  amrex::IntVect cell = bnd_bx.smallEnd();
295 
296  // Construct CC velocities
297  amrex::MultiFab mf_cc_vel;
298  auto ba = vars_new[ilev][Vars::cons].boxArray();
299  auto dm = vars_new[ilev][Vars::cons].DistributionMap();
300  mf_cc_vel.define(ba, dm, AMREX_SPACEDIM, amrex::IntVect(1,1,1));
301  average_face_to_cellcenter(mf_cc_vel,0,
302  amrex::Array<const amrex::MultiFab*,3>{&vars_new[ilev][Vars::xvel],
303  &vars_new[ilev][Vars::yvel],
304  &vars_new[ilev][Vars::zvel]});
305 
306  // Construct MultiFab holding requested variables
307  amrex::MultiFab mf_cc_data;
308  mf_cc_data.define(ba, dm, ncomp, 1);
309 
310  int mf_comp = 0;
311 
312  if (containerHasElement(m_varnames, "density")) {
313  amrex::MultiFab::Copy(mf_cc_data, vars_new[ilev][Vars::cons], Rho_comp, mf_comp, 1, 0);
314  mf_comp += 1;
315  }
316 
317  if (containerHasElement(m_varnames, "x_velocity")) {
318  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 0, mf_comp, 1, 0);
319  mf_comp += 1;
320  }
321  if (containerHasElement(m_varnames, "y_velocity")) {
322  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 1, mf_comp, 1, 0);
323  mf_comp += 1;
324  }
325  if (containerHasElement(m_varnames, "z_velocity")) {
326  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 2, mf_comp, 1, 0);
327  mf_comp += 1;
328  }
329 
330  if (containerHasElement(m_varnames, "magvel")) {
331 #ifdef _OPENMP
332 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
333 #endif
334  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
335  const amrex::Box& tbx = mfi.tilebox();
336  auto const& dfab = mf_cc_data.array(mfi);
337  auto const& vfab = mf_cc_vel.array(mfi);
338 
339  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
340  {
341  dfab(i,j,k,mf_comp) = std::sqrt(vfab(i,j,k,0)*vfab(i,j,k,0)
342  + vfab(i,j,k,1)*vfab(i,j,k,1)
343  + vfab(i,j,k,2)*vfab(i,j,k,2)) ;
344  });
345  }
346  mf_comp += 1;
347  }
348 
349  if (containerHasElement(m_varnames, "theta")) {
350 #ifdef _OPENMP
351 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
352 #endif
353  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
354  const amrex::Box& tbx = mfi.tilebox();
355  auto const& dfab = mf_cc_data.array(mfi);
356  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
357 
358  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
359  {
360  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoTheta_comp) / cfab(i,j,k,Rho_comp);
361  });
362  }
363  mf_comp += 1;
364  }
365 
366  if (containerHasElement(m_varnames, "qv")) {
367  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp,
368  "qv sampling requested but moisture components not present in state");
369  if (qv_comp >= 0) AMREX_ALWAYS_ASSERT(qv_comp == mf_comp);
370  qv_comp = mf_comp;
371 #ifdef _OPENMP
372 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
373 #endif
374  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
375  const amrex::Box& tbx = mfi.tilebox();
376  auto const& dfab = mf_cc_data.array(mfi);
377  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
378 
379  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
380  {
381  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
382  });
383  }
384  mf_comp += 1;
385  }
386  if (containerHasElement(m_varnames, "qc")) {
387  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ2_comp,
388  "qc sampling requested but moisture components not present in state");
389 #ifdef _OPENMP
390 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
391 #endif
392  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
393  const amrex::Box& tbx = mfi.tilebox();
394  auto const& dfab = mf_cc_data.array(mfi);
395  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
396 
397  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
398  {
399  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ2_comp) / cfab(i,j,k,Rho_comp);
400  });
401  }
402  mf_comp += 1;
403  }
404 
405  if (containerHasElement(m_varnames, "pressure")) {
406  if (vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp) {
407  // moist pressure: use qv from dfab if already sampled, else compute inline
408 #ifdef _OPENMP
409 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
410 #endif
411  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
412  const amrex::Box& tbx = mfi.tilebox();
413  auto const& dfab = mf_cc_data.array(mfi);
414  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
415 
416  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
417  {
418  amrex::Real qv_val = (qv_comp >= 0) ? dfab(i,j,k,qv_comp)
419  : cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
420  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp), qv_val);
421  });
422  }
423  } else {
424  // dry pressure
425 #ifdef _OPENMP
426 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
427 #endif
428  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
429  const amrex::Box& tbx = mfi.tilebox();
430  auto const& dfab = mf_cc_data.array(mfi);
431  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
432 
433  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
434  {
435  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp));
436  });
437  }
438  }
439  mf_comp += 1;
440  }
441 
442  m_lev[iline] = ilev;
443  m_ls_mf[iline] = get_line_data(mf_cc_data, dir, cell, bnd_bx);
444 
445  // We can stop if we got the entire line
446  auto min_bnd_bx = m_ls_mf[iline].boxArray().minimalBox();
447  if (bnd_bx == min_bnd_bx) { break; }
448 
449  } // ilev
450 
451  // Fill bnd_bx if use_rbx
452  if (m_use_real_bx) {
453  m_bnd_bx[iline] = m_ls_mf[iline].boxArray().minimalBox();
454  }
455 
456  }// iline
457  }
458 
459  void
460  write_sample_data (amrex::Vector<amrex::Real>& time,
461  amrex::Vector<int>& level_steps,
462  amrex::Vector<amrex::IntVect>& ref_ratio,
463  amrex::Vector<amrex::Geometry>& geom)
464  {
465  if (m_write_ascii) {
466  write_line_ascii(time);
467  } else {
468  write_line_plotfile(time, level_steps, ref_ratio, geom);
469  }
470  }
471 
472  void
473  write_line_ascii (amrex::Vector<amrex::Real>& time)
474  {
475  // same as definitions in ERF.H
476  constexpr int datwidth = 14;
477  constexpr int datprecision = 6;
478  constexpr int timeprecision = 13;
479 
480  int nline = static_cast<int>(m_ls_mf.size());
481  int nvar = static_cast<int>(m_varnames.size());
482  for (int iline(0); iline<nline; ++iline) {
483  int dir = m_dir[iline];
484  int lev = m_lev[iline];
485  amrex::Real m_time = time[lev];
486  amrex::Box m_dom = m_bnd_bx[iline];
487 
488  for (int ivar(0); ivar<nvar; ++ivar) {
489  // Convert multifab to vector
490  amrex::Gpu::HostVector<amrex::Real> vec = sumToLine(m_ls_mf[iline], ivar, 1, m_dom, dir);
491 
492  if (amrex::ParallelDescriptor::IOProcessor()) {
493  int ifile = iline*nvar + ivar;
494  std::ostream& fs = *m_datastream[ifile];
495  fs << std::setw(datwidth) << std::setprecision(timeprecision) << m_time
496  << std::setw(datwidth) << std::setprecision(datprecision);
497  for (const auto& val : vec) {
498  fs << " " << val;
499  }
500  fs << std::endl;
501  }
502  }
503  }
504  }
505 
506  void
507  write_line_plotfile (amrex::Vector<amrex::Real>& time,
508  amrex::Vector<int>& level_steps,
509  amrex::Vector<amrex::IntVect>& ref_ratio,
510  amrex::Vector<amrex::Geometry>& geom)
511  {
512  int nline = static_cast<int>(m_ls_mf.size());
513  for (int iline(0); iline<nline; ++iline) {
514  // Data members that can be used as-is
515  int dir = m_dir[iline];
516  int lev = m_lev[iline];
517  amrex::Real m_time = time[lev];
518  amrex::Vector<int> m_level_steps = {level_steps[lev]};
519  amrex::Vector<amrex::IntVect> m_ref_ratio = {ref_ratio[lev]};
520 
521  // Create modified geometry object corresponding to the line
522  auto plo = geom[lev].ProbLo();
523  auto dx = geom[lev].CellSize();
524  amrex::Vector<amrex::Geometry> m_geom; m_geom.resize(1);
525  amrex::Vector<int> is_per(AMREX_SPACEDIM,0);
526  amrex::Box m_dom = m_bnd_bx[iline];
527  amrex::RealBox m_rb;
528  for (int d(0); d<AMREX_SPACEDIM; ++d) {
529  amrex::Real offset = (d==dir) ? 0 : myhalf;
530  amrex::Real lo = plo[d] + ( m_dom.smallEnd(d) - offset ) * dx[d];
531  amrex::Real hi = plo[d] + ( m_dom.bigEnd(d) + offset ) * dx[d];
532 
533  m_rb.setLo(d,lo);
534  m_rb.setHi(d,hi);
535 
536  is_per[d] = geom[lev].isPeriodic(d);
537  }
538  m_geom[0].define(m_dom, &m_rb, geom[lev].Coord(), is_per.data());
539 
540  // Create plotfile name
541  std::string name_line = m_name[iline];
542  name_line += "_step_";
543  std::string plotfilename = amrex::Concatenate(name_line, m_level_steps[0], 5);
544 
545  // Get the data
546  amrex::Vector<const amrex::MultiFab*> mf = {&(m_ls_mf[iline])};
547 
548  // Write each line
549  WriteMultiLevelPlotfile(plotfilename, 1, mf,
550  m_varnames, m_geom, m_time,
551  m_level_steps, m_ref_ratio);
552  }
553  }
554 
555  amrex::Vector<int> m_dir;
556  amrex::Vector<int> m_lev;
557  amrex::Vector<amrex::Box> m_bnd_bx;
558  amrex::Vector<amrex::RealBox> m_bnd_rbx;
559  amrex::Vector<amrex::MultiFab> m_ls_mf;
560  amrex::Vector<std::string> m_name;
561 
562  bool m_use_real_bx{false};
563  bool m_write_ascii{false};
564  amrex::Vector<std::string> m_varnames {"magvel","theta"};
565  amrex::Vector<std::unique_ptr<std::fstream> > m_datastream;
566 };
567 
568 
570 {
572  {
573  amrex::ParmParse pp("erf");
574 
575  // Count number of lo and hi points define the plane
576  int n_plane_lo = pp.countval("sample_plane_lo") / AMREX_SPACEDIM;
577  int n_plane_hi = pp.countval("sample_plane_hi") / AMREX_SPACEDIM;
578  int n_plane_dir = pp.countval("sample_plane_dir");
579  AMREX_ALWAYS_ASSERT( (n_plane_lo==n_plane_hi ) &&
580  (n_plane_lo==n_plane_dir) );
581 
582  // Parse the data
583  if (n_plane_lo > 0) {
584  // Parse lo
585  amrex::Vector<amrex::Real> r_lo; r_lo.resize(n_plane_lo*AMREX_SPACEDIM);
586  amrex::Vector<amrex::Vector<amrex::Real>> rv_lo;
587  pp.queryarr("sample_plane_lo",r_lo,0,n_plane_lo*AMREX_SPACEDIM);
588  for (int i(0); i < n_plane_lo; i++) {
589  amrex::Vector<amrex::Real> rv = {r_lo[AMREX_SPACEDIM*i+0],
590  r_lo[AMREX_SPACEDIM*i+1],
591  r_lo[AMREX_SPACEDIM*i+2]};
592  rv_lo.push_back(rv);
593  }
594 
595  // Parse hi
596  amrex::Vector<amrex::Real> r_hi; r_hi.resize(n_plane_hi*AMREX_SPACEDIM);
597  amrex::Vector<amrex::Vector<amrex::Real>> rv_hi;
598  pp.queryarr("sample_plane_hi",r_hi,0,n_plane_hi*AMREX_SPACEDIM);
599  for (int i(0); i < n_plane_hi; i++) {
600  amrex::Vector<amrex::Real> rv = {r_hi[AMREX_SPACEDIM*i+0],
601  r_hi[AMREX_SPACEDIM*i+1],
602  r_hi[AMREX_SPACEDIM*i+2]};
603  rv_hi.push_back(rv);
604  }
605 
606  // Construct vector of bounding real boxes
607  m_bnd_rbx.resize(n_plane_lo);
608  for (int i(0); i < n_plane_hi; i++){
609  amrex::RealBox rbx(rv_lo[i].data(),rv_hi[i].data());
610  m_bnd_rbx[i] = rbx;
611  }
612 
613  // Parse directionality
614  m_dir.resize(n_plane_dir);
615  pp.queryarr("sample_plane_dir",m_dir,0,n_plane_dir);
616 
617  // Parse names
618  std::string name_base = "plt_plane_";
619  m_name.resize(n_plane_lo);
620  int n_names = pp.countval("sample_plane_name");
621  if (n_names > 0) {
622  AMREX_ALWAYS_ASSERT( n_names==n_plane_lo );
623  pp.queryarr("sample_plane_name",m_name,0,n_names);
624  } else {
625  for (int iplane(0); iplane<n_plane_lo; ++iplane) {
626  m_name[iplane] = amrex::Concatenate(name_base, iplane , 5);
627  }
628  }
629 
630  // Allocate space for level indicator
631  m_lev.resize(n_plane_dir,0);
632 
633  // Allocate space for MF pointers
634  m_ps_mf.resize(n_plane_lo);
635 
636  // Get requested vars
637  if (pp.countval("plane_sampling_vars") > 0) {
638  m_varnames.clear();
639  amrex::Vector<std::string> requested_vars;
640  pp.queryarr("plane_sampling_vars",requested_vars);
641  amrex::Print() << "Selected plane sampling vars :";
642  if (containerHasElement(requested_vars, "density")) {
643  m_varnames.push_back("density");
644  amrex::Print() << " " << "density";
645  }
646  if (containerHasElement(requested_vars, "x_velocity")) {
647  m_varnames.push_back("x_velocity");
648  amrex::Print() << " " << "x_velocity";
649  }
650  if (containerHasElement(requested_vars, "y_velocity")) {
651  m_varnames.push_back("y_velocity");
652  amrex::Print() << " " << "y_velocity";
653  }
654  if (containerHasElement(requested_vars, "z_velocity")) {
655  m_varnames.push_back("z_velocity");
656  amrex::Print() << " " << "z_velocity";
657  }
658  if (containerHasElement(requested_vars, "magvel")) {
659  m_varnames.push_back("magvel");
660  amrex::Print() << " " << "magvel";
661  }
662  if (containerHasElement(requested_vars, "theta")) {
663  m_varnames.push_back("theta");
664  amrex::Print() << " " << "theta";
665  }
666  if (containerHasElement(requested_vars, "qv")) {
667  m_varnames.push_back("qv");
668  amrex::Print() << " " << "qv";
669  }
670  if (containerHasElement(requested_vars, "qc")) {
671  m_varnames.push_back("qc");
672  amrex::Print() << " " << "qc";
673  }
674  if (containerHasElement(requested_vars, "pressure")) {
675  m_varnames.push_back("pressure");
676  amrex::Print() << " " << "pressure";
677  }
678  amrex::Print() << std::endl;
679  }
680  }
681  }
682 
683  // This must match what is in AMReX_MultiFabUtil.H
684  amrex::Box
685  getIndexBox (const amrex::RealBox& real_box,
686  const amrex::Geometry& geom) {
687  amrex::IntVect slice_lo, slice_hi;
688 
689  AMREX_D_TERM(slice_lo[0]=static_cast<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
690  slice_lo[1]=static_cast<int>(std::floor((real_box.lo(1) - geom.ProbLo(1))/geom.CellSize(1)));,
691  slice_lo[2]=static_cast<int>(std::floor((real_box.lo(2) - geom.ProbLo(2))/geom.CellSize(2))););
692 
693  AMREX_D_TERM(slice_hi[0]=static_cast<int>(std::floor((real_box.hi(0) - geom.ProbLo(0))/geom.CellSize(0)));,
694  slice_hi[1]=static_cast<int>(std::floor((real_box.hi(1) - geom.ProbLo(1))/geom.CellSize(1)));,
695  slice_hi[2]=static_cast<int>(std::floor((real_box.hi(2) - geom.ProbLo(2))/geom.CellSize(2))););
696 
697  return amrex::Box(slice_lo, slice_hi) & geom.Domain();
698  }
699 
700  void
701  get_sample_data (amrex::Vector<amrex::Geometry>& geom,
702  amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new)
703  {
704  int nlev = static_cast<int>(vars_new.size());
705  int nplane = static_cast<int>(m_bnd_rbx.size());
706  int ncomp = static_cast<int>(m_varnames.size());
707  bool interpolate = true;
708 
709  int qv_comp = -1;
710 
711  // Loop over each plane
712  for (int iplane(0); iplane<nplane; ++iplane) {
713  int dir = m_dir[iplane];
714  amrex::RealBox bnd_rbx = m_bnd_rbx[iplane];
715  amrex::Real point = bnd_rbx.lo(dir);
716 
717  // Search each level to get the finest data possible
718  for (int ilev(nlev-1); ilev>=0; --ilev) {
719 
720  // Construct CC velocities
721  amrex::MultiFab mf_cc_vel;
722  auto ba = vars_new[ilev][Vars::cons].boxArray();
723  auto dm = vars_new[ilev][Vars::cons].DistributionMap();
724  mf_cc_vel.define(ba, dm, AMREX_SPACEDIM, amrex::IntVect(1,1,1));
725  average_face_to_cellcenter(mf_cc_vel,0,
726  amrex::Array<const amrex::MultiFab*,3>{&vars_new[ilev][Vars::xvel],
727  &vars_new[ilev][Vars::yvel],
728  &vars_new[ilev][Vars::zvel]});
729 
730  // Construct MultiFab holding requested variables
731  amrex::MultiFab mf_cc_data;
732  mf_cc_data.define(ba, dm, ncomp, 1);
733 
734  int mf_comp = 0;
735 
736  if (containerHasElement(m_varnames, "density")) {
737  amrex::MultiFab::Copy(mf_cc_data, vars_new[ilev][Vars::cons], Rho_comp, mf_comp, 1, 0);
738  mf_comp += 1;
739  }
740 
741  if (containerHasElement(m_varnames, "x_velocity")) {
742  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 0, mf_comp, 1, 0);
743  mf_comp += 1;
744  }
745  if (containerHasElement(m_varnames, "y_velocity")) {
746  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 1, mf_comp, 1, 0);
747  mf_comp += 1;
748  }
749  if (containerHasElement(m_varnames, "z_velocity")) {
750  amrex::MultiFab::Copy(mf_cc_data, mf_cc_vel, 2, mf_comp, 1, 0);
751  mf_comp += 1;
752  }
753 
754  if (containerHasElement(m_varnames, "magvel")) {
755 #ifdef _OPENMP
756 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
757 #endif
758  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
759  const amrex::Box& tbx = mfi.tilebox();
760  auto const& dfab = mf_cc_data.array(mfi);
761  auto const& vfab = mf_cc_vel.array(mfi);
762 
763  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
764  {
765  dfab(i,j,k,mf_comp) = std::sqrt(vfab(i,j,k,0)*vfab(i,j,k,0)
766  + vfab(i,j,k,1)*vfab(i,j,k,1)
767  + vfab(i,j,k,2)*vfab(i,j,k,2));
768  });
769  }
770  mf_comp += 1;
771  }
772 
773  if (containerHasElement(m_varnames, "theta")) {
774 #ifdef _OPENMP
775 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
776 #endif
777  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
778  const amrex::Box& tbx = mfi.tilebox();
779  auto const& dfab = mf_cc_data.array(mfi);
780  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
781 
782  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
783  {
784  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoTheta_comp) / cfab(i,j,k,Rho_comp);
785  });
786  }
787  mf_comp += 1;
788  }
789 
790  if (containerHasElement(m_varnames, "qv")) {
791  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp,
792  "qv sampling requested but moisture components not present in state");
793  if (qv_comp >= 0) AMREX_ALWAYS_ASSERT(qv_comp == mf_comp);
794  qv_comp = mf_comp;
795 #ifdef _OPENMP
796 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
797 #endif
798  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
799  const amrex::Box& tbx = mfi.tilebox();
800  auto const& dfab = mf_cc_data.array(mfi);
801  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
802 
803  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
804  {
805  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
806  });
807  }
808  mf_comp += 1;
809  }
810  if (containerHasElement(m_varnames, "qc")) {
811  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(vars_new[ilev][Vars::cons].nComp() > RhoQ2_comp,
812  "qc sampling requested but moisture components not present in state");
813 #ifdef _OPENMP
814 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
815 #endif
816  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
817  const amrex::Box& tbx = mfi.tilebox();
818  auto const& dfab = mf_cc_data.array(mfi);
819  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
820 
821  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
822  {
823  dfab(i,j,k,mf_comp) = cfab(i,j,k,RhoQ2_comp) / cfab(i,j,k,Rho_comp);
824  });
825  }
826  mf_comp += 1;
827  }
828 
829  if (containerHasElement(m_varnames, "pressure")) {
830  if (vars_new[ilev][Vars::cons].nComp() > RhoQ1_comp) {
831  // moist pressure: use qv from dfab if already sampled, else compute inline
832 #ifdef _OPENMP
833 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
834 #endif
835  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
836  const amrex::Box& tbx = mfi.tilebox();
837  auto const& dfab = mf_cc_data.array(mfi);
838  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
839 
840  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
841  {
842  amrex::Real qv_val = (qv_comp >= 0) ? dfab(i,j,k,qv_comp)
843  : cfab(i,j,k,RhoQ1_comp) / cfab(i,j,k,Rho_comp);
844  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp), qv_val);
845  });
846  }
847  } else {
848  // dry pressure
849 #ifdef _OPENMP
850 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
851 #endif
852  for (amrex::MFIter mfi(mf_cc_data, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
853  const amrex::Box& tbx = mfi.tilebox();
854  auto const& dfab = mf_cc_data.array(mfi);
855  auto const& cfab = vars_new[ilev][Vars::cons].array(mfi);
856 
857  amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
858  {
859  dfab(i,j,k,mf_comp) = getPgivenRTh(cfab(i,j,k,RhoTheta_comp));
860  });
861  }
862  }
863  mf_comp += 1;
864  }
865 
866  m_lev[iplane] = ilev;
867  m_ps_mf[iplane] = get_slice_data(dir, point, mf_cc_data, geom[ilev],
868  0, ncomp, interpolate, bnd_rbx);
869 
870  // We can stop if we got the entire plane
871  auto min_bnd_bx = m_ps_mf[iplane]->boxArray().minimalBox();
872  amrex::Box bnd_bx = getIndexBox(bnd_rbx, geom[ilev]);
873  if (bnd_bx == min_bnd_bx) { break; }
874 
875  } // ilev
876  }// iplane
877  }
878 
879  void
880  write_sample_data (amrex::Vector<amrex::Real>& time,
881  amrex::Vector<int>& level_steps,
882  amrex::Vector<amrex::IntVect>& ref_ratio,
883  amrex::Vector<amrex::Geometry>& geom)
884  {
885  int nplane = m_ps_mf.size();
886  for (int iplane(0); iplane<nplane; ++iplane) {
887  // Data members that can be used as-is
888  int dir = m_dir[iplane];
889  int lev = m_lev[iplane];
890  amrex::Real m_time = time[lev];
891  amrex::Vector<int> m_level_steps = {level_steps[lev]};
892  amrex::Vector<amrex::IntVect> m_ref_ratio = {ref_ratio[lev]};
893 
894  // Create modified geometry object corresponding to the plane
895  amrex::RealBox m_rb = m_bnd_rbx[iplane];
896  amrex::Box m_dom = getIndexBox(m_rb, geom[lev]);
897  amrex::Real point = m_rb.hi(dir);
898  amrex::Vector<int> is_per(AMREX_SPACEDIM,0);
899  for (int d(0); d<AMREX_SPACEDIM; ++d) {
900  if (d==dir) {
901  m_rb.setLo(d,point-myhalf*geom[lev].CellSize(d));
902  m_rb.setHi(d,point+myhalf*geom[lev].CellSize(d));
903  }
904  is_per[d] = geom[lev].isPeriodic(d);
905  }
906  amrex::Vector<amrex::Geometry> m_geom; m_geom.resize(1);
907  m_geom[0].define(m_dom, &m_rb, geom[lev].Coord(), is_per.data());
908 
909  // Create plotfile name
910  std::string name_plane = m_name[iplane];
911  name_plane += "_step_";
912  std::string plotfilename = amrex::Concatenate(name_plane, m_level_steps[0], 5);
913 
914  // Get the data
915  amrex::Vector<const amrex::MultiFab*> mf = {m_ps_mf[iplane].get()};
916 
917  // Write each plane
918  WriteMultiLevelPlotfile(plotfilename, 1, mf,
919  m_varnames, m_geom, m_time,
920  m_level_steps, m_ref_ratio);
921  } // iplane
922  }
923 
924  amrex::Vector<int> m_dir;
925  amrex::Vector<int> m_lev;
926  amrex::Vector<amrex::RealBox> m_bnd_rbx;
927  amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_ps_mf;
928  amrex::Vector<std::string> m_name;
929 
930  amrex::Vector<std::string> m_varnames {"magvel","theta"};
931 };
932 #endif
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
bool containerHasElement(const V &iterable, const T &query)
Definition: ERF_Container.H:5
Coord
Definition: ERF_DataStruct.H:92
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
struct @28 out
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
ParmParse pp("prob")
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor(grown_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);})
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
amrex::Real Real
Definition: ERF_ShocInterface.H:19
@ fs
Definition: ERF_AdvanceMorrison.cpp:122
@ xvel
Definition: ERF_IndexDefines.H:175
@ cons
Definition: ERF_IndexDefines.H:174
@ zvel
Definition: ERF_IndexDefines.H:177
@ yvel
Definition: ERF_IndexDefines.H:176
Definition: ERF_SampleData.H:18
amrex::Vector< std::string > m_varnames
Definition: ERF_SampleData.H:564
amrex::Vector< amrex::MultiFab > m_ls_mf
Definition: ERF_SampleData.H:559
void write_coords(amrex::Vector< std::unique_ptr< amrex::MultiFab > > &z_phys_cc, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:223
amrex::Vector< std::unique_ptr< std::fstream > > m_datastream
Definition: ERF_SampleData.H:565
void write_line_plotfile(amrex::Vector< amrex::Real > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:507
bool m_use_real_bx
Definition: ERF_SampleData.H:562
amrex::Vector< int > m_dir
Definition: ERF_SampleData.H:555
amrex::Vector< std::string > m_name
Definition: ERF_SampleData.H:560
void write_line_ascii(amrex::Vector< amrex::Real > &time)
Definition: ERF_SampleData.H:473
LineSampler()
Definition: ERF_SampleData.H:19
amrex::Vector< amrex::Box > m_bnd_bx
Definition: ERF_SampleData.H:557
void get_sample_data(amrex::Vector< amrex::Geometry > &geom, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_new)
Definition: ERF_SampleData.H:276
void write_sample_data(amrex::Vector< amrex::Real > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:460
amrex::Box getIndexBox(const amrex::RealBox &real_box, const amrex::Geometry &geom)
Definition: ERF_SampleData.H:207
bool m_write_ascii
Definition: ERF_SampleData.H:563
amrex::Vector< amrex::RealBox > m_bnd_rbx
Definition: ERF_SampleData.H:558
amrex::Vector< int > m_lev
Definition: ERF_SampleData.H:556
Definition: ERF_SampleData.H:570
void write_sample_data(amrex::Vector< amrex::Real > &time, amrex::Vector< int > &level_steps, amrex::Vector< amrex::IntVect > &ref_ratio, amrex::Vector< amrex::Geometry > &geom)
Definition: ERF_SampleData.H:880
amrex::Vector< int > m_dir
Definition: ERF_SampleData.H:924
amrex::Box getIndexBox(const amrex::RealBox &real_box, const amrex::Geometry &geom)
Definition: ERF_SampleData.H:685
amrex::Vector< int > m_lev
Definition: ERF_SampleData.H:925
amrex::Vector< std::string > m_name
Definition: ERF_SampleData.H:928
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_ps_mf
Definition: ERF_SampleData.H:927
amrex::Vector< amrex::RealBox > m_bnd_rbx
Definition: ERF_SampleData.H:926
PlaneSampler()
Definition: ERF_SampleData.H:571
amrex::Vector< std::string > m_varnames
Definition: ERF_SampleData.H:930
void get_sample_data(amrex::Vector< amrex::Geometry > &geom, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_new)
Definition: ERF_SampleData.H:701