REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_Plotfile.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2#include "AMReX_Interp_3D_C.H"
3#include "AMReX_PlotFileUtil.H"
4
5using namespace amrex;
6
7PhysBCFunctNoOp null_bc_for_fill;
8
9template<typename V, typename T>
10bool containerHasElement(const V& iterable, const T& query) {
11 return std::find(iterable.begin(), iterable.end(), query) != iterable.end();
12}
13
14/**
15 * @param pp_plot_var_names list of variable names to plot read in from parameter file
16 */
17void
18REMORA::setPlotVariables (const std::string& pp_plot_var_names)
19{
20 ParmParse pp(pp_prefix);
21
22 if (pp.contains(pp_plot_var_names.c_str()))
23 {
24 std::string nm;
25
26 int nPltVars = pp.countval(pp_plot_var_names.c_str());
27
28 for (int i = 0; i < nPltVars; i++)
29 {
30 pp.get(pp_plot_var_names.c_str(), nm, i);
31
32 // Add the named variable to our list of plot variables
33 // if it is not already in the list
34 if (!containerHasElement(plot_var_names, nm)) {
35 plot_var_names.push_back(nm);
36 }
37 }
38 } else {
39 //
40 // The default is to add none of the variables to the list
41 //
42 plot_var_names.clear();
43 }
44
45 // Get state variables in the same order as we define them,
46 // since they may be in any order in the input list
47 Vector<std::string> tmp_plot_names;
48
49 for (int i = 0; i < NCONS; ++i) {
50 if ( containerHasElement(plot_var_names, cons_names[i]) ) {
51 tmp_plot_names.push_back(cons_names[i]);
52 }
53 }
54 // Check for velocity since it's not in cons_names
55 // If we are asked for any velocity component, we will need them all
56 if (containerHasElement(plot_var_names, "x_velocity") ||
57 containerHasElement(plot_var_names, "y_velocity") ||
58 containerHasElement(plot_var_names, "z_velocity")) {
59 tmp_plot_names.push_back("x_velocity");
60 tmp_plot_names.push_back("y_velocity");
61 tmp_plot_names.push_back("z_velocity");
62 }
63
64 // If we are asked for any location component, we will provide them all
65 if (containerHasElement(plot_var_names, "x_cc") ||
66 containerHasElement(plot_var_names, "y_cc") ||
67 containerHasElement(plot_var_names, "z_cc")) {
68 tmp_plot_names.push_back("x_cc");
69 tmp_plot_names.push_back("y_cc");
70 tmp_plot_names.push_back("z_cc");
71 }
72
73 for (int i = 0; i < derived_names.size(); ++i) {
74 if ( containerHasElement(plot_var_names, derived_names[i]) ) {
75 tmp_plot_names.push_back(derived_names[i]);
76 } // if
77 } // i
78
79#ifdef REMORA_USE_PARTICLES
80 const auto& particles_namelist( particleData.getNamesUnalloc() );
81 for (auto it = particles_namelist.cbegin(); it != particles_namelist.cend(); ++it) {
82 std::string tmp( (*it)+"_count" );
83 if (containerHasElement(plot_var_names, tmp) ) {
84 tmp_plot_names.push_back(tmp);
85 }
86 }
87#endif
88
89 // Check to see if we found all the requested variables
90 for (auto plot_name : plot_var_names) {
91 if (!containerHasElement(tmp_plot_names, plot_name)) {
92 Warning("\nWARNING: Requested to plot variable '" + plot_name + "' but it is not available");
93 }
94 }
95 plot_var_names = tmp_plot_names;
96}
97
98/**
99 * @param pp_plot_var_names variables to add to plot list
100 */
101void
102REMORA::appendPlotVariables (const std::string& pp_plot_var_names)
103{
104 ParmParse pp(pp_prefix);
105
106 if (pp.contains(pp_plot_var_names.c_str())) {
107 std::string nm;
108 int nPltVars = pp.countval(pp_plot_var_names.c_str());
109 for (int i = 0; i < nPltVars; i++) {
110 pp.get(pp_plot_var_names.c_str(), nm, i);
111 // Add the named variable to our list of plot variables
112 // if it is not already in the list
113 if (!containerHasElement(plot_var_names, nm)) {
114 plot_var_names.push_back(nm);
115 }
116 }
117 }
118
119 Vector<std::string> tmp_plot_names(0);
120#ifdef REMORA_USE_PARTICLES
121 Vector<std::string> particle_mesh_plot_names;
122 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
123 for (int i = 0; i < particle_mesh_plot_names.size(); i++) {
124 std::string tmp(particle_mesh_plot_names[i]);
125 if (containerHasElement(plot_var_names, tmp) ) {
126 tmp_plot_names.push_back(tmp);
127 }
128 }
129#endif
130
131 for (int i = 0; i < tmp_plot_names.size(); i++) {
132 plot_var_names.push_back( tmp_plot_names[i] );
133 }
134
135 // Finally, check to see if we found all the requested variables
136 for (const auto& plot_name : plot_var_names) {
137 if (!containerHasElement(plot_var_names, plot_name)) {
138 if (amrex::ParallelDescriptor::IOProcessor()) {
139 Warning("\nWARNING: Requested to plot variable '" + plot_name + "' but it is not available");
140 }
141 }
142 }
143}
144
145// Write plotfile to disk
146void
147REMORA::WritePlotFile ()
148{
149 Vector<std::string> varnames;
150 varnames.insert(varnames.end(), plot_var_names.begin(), plot_var_names.end());
151
152 const int ncomp_mf = varnames.size();
153 const auto ngrow_vars = IntVect(NGROW-1,NGROW-1,0);
154
155 if (ncomp_mf == 0) {
156 return;
157 }
158
159 // We fillpatch here because some of the derived quantities require derivatives
160 // which require ghost cells to be filled. Don't fill the boundary, though.
161 for (int lev = 0; lev <= finest_level; ++lev) {
162 FillPatchNoBC(lev, t_new[lev], *cons_new[lev], cons_new, BdyVars::t,0,true,false);
163 FillPatchNoBC(lev, t_new[lev], *xvel_new[lev], xvel_new, BdyVars::u,0,true,false);
164 FillPatchNoBC(lev, t_new[lev], *yvel_new[lev], yvel_new, BdyVars::v,0,true,false);
165 FillPatchNoBC(lev, t_new[lev], *zvel_new[lev], zvel_new, BdyVars::null,0,true,false);
166 }
167
168 Real fill_value = 0.0_rt;
169 for (int lev = 0; lev <= finest_level; ++lev) {
170 mask_arrays_for_write(lev, (Real) fill_value, 0.0_rt);
171 }
172
173 // Array of MultiFabs to hold the plotfile data
174 Vector<MultiFab> mf(finest_level+1);
175 for (int lev = 0; lev <= finest_level; ++lev) {
176 mf[lev].define(grids[lev], dmap[lev], ncomp_mf, ngrow_vars);
177 }
178
179 // Array of MultiFabs for nodal data
180 Vector<MultiFab> mf_nd(finest_level+1);
181 for (int lev = 0; lev <= finest_level; ++lev) {
182 BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes();
183 mf_nd[lev].define(nodal_grids, dmap[lev], AMREX_SPACEDIM, 0);
184 mf_nd[lev].setVal(0.);
185 }
186
187 // Vector of MultiFabs for face-centered velocity
188 Vector<MultiFab> mf_u(finest_level+1);
189 Vector<MultiFab> mf_v(finest_level+1);
190 Vector<MultiFab> mf_w(finest_level+1);
191 if (plot_staggered_vels) {
192 for (int lev = 0; lev <= finest_level; ++lev) {
193 BoxArray grid_stag_u(grids[lev]); grid_stag_u.surroundingNodes(0);
194 BoxArray grid_stag_v(grids[lev]); grid_stag_v.surroundingNodes(1);
195 BoxArray grid_stag_w(grids[lev]); grid_stag_w.surroundingNodes(2);
196 mf_u[lev].define(grid_stag_u, dmap[lev], 1, 0);
197 mf_v[lev].define(grid_stag_v, dmap[lev], 1, 0);
198 mf_w[lev].define(grid_stag_w, dmap[lev], 1, 0);
199 MultiFab::Copy(mf_u[lev],*xvel_new[lev],0,0,1,0);
200 MultiFab::Copy(mf_v[lev],*yvel_new[lev],0,0,1,0);
201 MultiFab::Copy(mf_w[lev],*zvel_new[lev],0,0,1,0);
202 }
203 }
204
205 // Array of MultiFabs for cell-centered velocity
206 Vector<MultiFab> mf_cc_vel(finest_level+1);
207
208 if (containerHasElement(plot_var_names, "x_velocity") ||
209 containerHasElement(plot_var_names, "y_velocity") ||
210 containerHasElement(plot_var_names, "z_velocity") ||
211 containerHasElement(plot_var_names, "vorticity") ) {
212
213 for (int lev = 0; lev <= finest_level; ++lev) {
214 mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,0));
215 mf_cc_vel[lev].setVal(0.0_rt); // zero out velocity in case we have any wall boundaries
216 average_face_to_cellcenter(mf_cc_vel[lev],0,
217 Array<const MultiFab*,3>{xvel_new[lev],yvel_new[lev],zvel_new[lev]},IntVect(1,1,0));
218 mf_cc_vel[lev].FillBoundary(geom[lev].periodicity());
219 } // lev
220
221 // We need ghost cells if computing vorticity
222 amrex::Interpolater* mapper = &cell_cons_interp;
223 if ( containerHasElement(plot_var_names, "vorticity") ) {
224 for (int lev = 1; lev <= finest_level; ++lev) {
225 Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
226 Vector<Real> ftime = {t_new[lev], t_new[lev]};
227 Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
228 Vector<Real> ctime = {t_new[lev], t_new[lev]};
229
230 MultiFab mf_to_fill;
231 amrex::FillPatchTwoLevels(mf_cc_vel[lev], mf_cc_vel[lev].nGrowVect(), IntVect(0,0,0),
232 t_new[lev], cmf, ctime, fmf, ftime,
233 0, 0, mf_cc_vel[lev].nComp(), geom[lev-1], geom[lev],
234 refRatio(lev-1), mapper, domain_bcs_type, BCVars::foextrap_bc);
235 } // lev
236 } // if
237 } // if
238
239 for (int lev = 0; lev <= finest_level; ++lev)
240 {
241 int mf_comp = 0;
242
243 // First, copy any of the conserved state variables into the output plotfile
244 AMREX_ALWAYS_ASSERT(cons_names.size() == NCONS);
245 for (int i = 0; i < NCONS; ++i) {
246 if (containerHasElement(plot_var_names, cons_names[i])) {
247 if (cons_new[lev]->contains_nan() || cons_new[lev]->contains_inf()) {
248 amrex::Abort("Found while writing output: Cons (salt, temp, or scalar, etc) contains nan or inf");
249 }
250 MultiFab::Copy(mf[lev],*cons_new[lev],i,mf_comp,1,ngrow_vars);
251 mf_comp++;
252 }
253 } // NCONS
254
255 // Next, check for velocities
256 if (containerHasElement(plot_var_names, "x_velocity")) {
257 if (mf_cc_vel[lev].contains_nan(0,1) || mf_cc_vel[lev].contains_inf(0,1)) {
258 amrex::Abort("Found while writing output: u velocity contains nan or inf");
259 }
260 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 0, mf_comp, 1, 0);
261 mf_comp += 1;
262 }
263 if (containerHasElement(plot_var_names, "y_velocity")) {
264 if (mf_cc_vel[lev].contains_nan(1,1) || mf_cc_vel[lev].contains_inf(1,1)) {
265 amrex::Abort("Found while writing output: v velocity contains nan or inf");
266 }
267 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 1, mf_comp, 1, 0);
268 mf_comp += 1;
269 }
270 if (containerHasElement(plot_var_names, "z_velocity")) {
271 if (mf_cc_vel[lev].contains_nan(2,1) || mf_cc_vel[lev].contains_inf(2,1)) {
272 amrex::Abort("Found while writing output: z velocity contains nan or inf");
273 }
274 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 2, mf_comp, 1, 0);
275 mf_comp += 1;
276 }
277
278 // Define standard process for calling the functions in Derive.cpp
279 auto calculate_derived = [&](const std::string& der_name,
280 decltype(derived::remora_dernull)& der_function)
281 {
282 if (containerHasElement(plot_var_names, der_name)) {
283 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
284#ifdef _OPENMP
285#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
286#endif
287 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
288 {
289 const Box& bx = mfi.tilebox();
290 auto& dfab = dmf[mfi];
291
292 if (der_name == "vorticity") {
293 auto const& sfab = mf_cc_vel[lev][mfi];
294 der_function(bx, dfab, 0, 1, sfab, vec_pm[lev]->const_array(mfi), vec_pn[lev]->const_array(mfi), Geom(lev), t_new[0], nullptr, lev);
295 } else {
296 auto const& sfab = (*cons_new[lev])[mfi];
297 der_function(bx, dfab, 0, 1, sfab, vec_pm[lev]->const_array(mfi), vec_pn[lev]->const_array(mfi), Geom(lev), t_new[0], nullptr, lev);
298 }
299 }
300
301 mf_comp++;
302 }
303 };
304
305 // Note: All derived variables must be computed in order of "derived_names" defined in REMORA.H
306 calculate_derived("vorticity", derived::remora_dervort);
307
308 // Fill cell-centered location
309 Real dx = Geom()[lev].CellSizeArray()[0];
310 Real dy = Geom()[lev].CellSizeArray()[1];
311
312 // Next, check for location names -- if we write one we write all
313 if (containerHasElement(plot_var_names, "x_cc") ||
314 containerHasElement(plot_var_names, "y_cc") ||
315 containerHasElement(plot_var_names, "z_cc"))
316 {
317 MultiFab dmf(mf[lev], make_alias, mf_comp, AMREX_SPACEDIM);
318#ifdef _OPENMP
319#pragma omp parallel if (Gpu::notInLaunchRegion())
320#endif
321 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
322 const Box& bx = mfi.tilebox();
323 const Array4<Real> loc_arr = dmf.array(mfi);
324 const Array4<Real const> zp_arr = vec_z_phys_nd[lev]->const_array(mfi);
325
326 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
327 loc_arr(i,j,k,0) = (i+0.5_rt) * dx;
328 loc_arr(i,j,k,1) = (j+0.5_rt) * dy;
329 loc_arr(i,j,k,2) = 0.125_rt * (zp_arr(i,j ,k ) + zp_arr(i+1,j ,k ) +
330 zp_arr(i,j+1,k ) + zp_arr(i+1,j+1,k ) +
331 zp_arr(i,j ,k+1) + zp_arr(i+1,j ,k+1) +
332 zp_arr(i,j+1,k+1) + zp_arr(i+1,j+1,k+1) );
333 });
334 } // mfi
335 mf_comp += AMREX_SPACEDIM;
336 } // if containerHasElement
337
338#ifdef REMORA_USE_PARTICLES
339 const auto& particles_namelist( particleData.getNames() );
340 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
341 if (containerHasElement(plot_var_names, std::string(particles_namelist[i]+"_count"))) {
342 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
343 temp_dat.setVal(0);
344 particleData[particles_namelist[i]]->Increment(temp_dat, lev);
345 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
346 mf_comp += 1;
347 }
348 }
349
350 Vector<std::string> particle_mesh_plot_names(0);
351 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
352 for (int i = 0; i < particle_mesh_plot_names.size(); i++) {
353 std::string plot_var_name(particle_mesh_plot_names[i]);
354 if (containerHasElement(plot_var_names, plot_var_name) ) {
355 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
356 temp_dat.setVal(0);
357 particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
358 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
359 mf_comp += 1;
360 }
361 }
362#endif
363
364 MultiFab::Copy(mf_nd[lev],*vec_z_phys_nd[lev],0,2,1,0);
365 Real dz = Geom()[lev].CellSizeArray()[2];
366 int N = Geom()[lev].Domain().size()[2];
367
368#ifdef _OPENMP
369#pragma omp parallel if (Gpu::notInLaunchRegion())
370#endif
371 for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
372 {
373 const Box& bx = mfi.tilebox();
374 Array4<Real> mf_arr = mf_nd[lev].array(mfi);
375 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
376 mf_arr(i,j,k,2) = mf_arr(i,j,k,2) + (N-k) * dz;
377 });
378 } // mfi
379
380 } // lev
381
382 std::string plotfilename = Concatenate(plot_file_name, istep[0], file_min_digits);
383
384 if (finest_level == 0)
385 {
386 if (plotfile_type == PlotfileType::amrex) {
387 amrex::Print() << "Writing plotfile " << plotfilename << "\n";
388 WriteMultiLevelPlotfileWithBathymetry(plotfilename, finest_level+1,
389 GetVecOfConstPtrs(mf),
390 GetVecOfConstPtrs(mf_nd),
391 GetVecOfConstPtrs(mf_u),
392 GetVecOfConstPtrs(mf_v),
393 GetVecOfConstPtrs(mf_w),
394 varnames,
395 Geom(),
396 t_new[0], istep, refRatio());
397 writeJobInfo(plotfilename);
398
399#ifdef REMORA_USE_PARTICLES
400 particleData.Checkpoint(plotfilename);
401#endif
402
403#ifdef REMORA_USE_HDF5
404 } else if (plotfile_type == PlotfileType::hdf5) {
405 amrex::Print() << "Writing plotfile " << plotfilename+"d01.h5" << "\n";
406 WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1,
407 GetVecOfConstPtrs(mf),
408 varnames,
409 Geom(), t_new[0], istep, refRatio());
410#endif
411 } else if (!(plotfile_type == PlotfileType::netcdf)) {
412 amrex::Abort("User specified unknown plot_filetype");
413 }
414
415 } else { // multilevel
416 if (plotfile_type == PlotfileType::amrex) {
417 amrex::Print() << "Writing plotfile " << plotfilename << "\n";
418 int lev0 = 0;
419 [[maybe_unused]] int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]);
420 bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 1) || (ref_ratio[lev0][1] == 1) ) ||
421 (ref_ratio[lev0][2] == 1) );
422 for (int lev = 1; lev < finest_level; lev++) {
423 any_ratio_one = any_ratio_one ||
424 ( ( (ref_ratio[lev][0] == 1) || (ref_ratio[lev][1] == 1) ) ||
425 (ref_ratio[lev][2] == 1) );
426 }
427 if (any_ratio_one && expand_plotvars_to_unif_rr) {
428 Vector<IntVect> r2(finest_level);
429 Vector<Geometry> g2(finest_level+1);
430 Vector<MultiFab> mf2(finest_level+1);
431
432 mf2[0].define(grids[0], dmap[0], ncomp_mf, 0);
433
434 // Copy level 0 as is
435 MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0);
436
437 // Define a new multi-level array of Geometry's so that we pass the new "domain" at lev > 0
438 Array<int,AMREX_SPACEDIM> periodicity =
439 {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)};
440 g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data());
441
442 r2[0] = IntVect(1,1,ref_ratio[0][0]);
443 for (int lev = 1; lev <= finest_level; ++lev) {
444 if (lev > 1) {
445 r2[lev-1][0] = 1;
446 r2[lev-1][1] = 1;
447 r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0];
448 }
449
450 mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
451
452 // Set the new problem domain
453 Box d2(Geom()[lev].Domain());
454 d2.refine(r2[lev-1]);
455
456 g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data());
457 }
458
459 // Make a vector of BCRec with default values so we can use it here -- note the values
460 // aren't actually used because we do PCInterp
461 amrex::Vector<amrex::BCRec> null_dom_bcs;
462 null_dom_bcs.resize(mf2[0].nComp());
463 for (int n = 0; n < mf2[0].nComp(); n++) {
464 for (int dir = 0; dir < AMREX_SPACEDIM; dir++) {
465 null_dom_bcs[n].setLo(dir, REMORABCType::int_dir);
466 null_dom_bcs[n].setHi(dir, REMORABCType::int_dir);
467 }
468 }
469
470 // Do piecewise interpolation of mf into mf2
471 for (int lev = 1; lev <= finest_level; ++lev) {
472 Interpolater* mapper_c = &pc_interp;
473 InterpFromCoarseLevel(mf2[lev], t_new[lev], mf[lev],
474 0, 0, mf2[lev].nComp(),
475 geom[lev], g2[lev],
476 null_bc_for_fill, 0, null_bc_for_fill, 0,
477 r2[lev-1], mapper_c, null_dom_bcs, 0);
478 }
479
480 // Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile
481 Vector<IntVect> rr(finest_level);
482 for (int lev = 0; lev < finest_level; ++lev) {
483 rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]);
484 }
485
486 WriteMultiLevelPlotfileWithBathymetry(plotfilename, finest_level+1,
487 GetVecOfConstPtrs(mf2),
488 GetVecOfConstPtrs(mf_nd),
489 GetVecOfConstPtrs(mf_u),
490 GetVecOfConstPtrs(mf_v),
491 GetVecOfConstPtrs(mf_w),
492 varnames,
493 g2,
494 t_new[0], istep, rr);
495 writeJobInfo(plotfilename);
496
497#ifdef REMORA_USE_PARTICLES
498 particleData.Checkpoint(plotfilename);
499#endif
500 } else {
501 WriteMultiLevelPlotfileWithBathymetry(plotfilename, finest_level+1,
502 GetVecOfConstPtrs(mf),
503 GetVecOfConstPtrs(mf_nd),
504 GetVecOfConstPtrs(mf_u),
505 GetVecOfConstPtrs(mf_v),
506 GetVecOfConstPtrs(mf_w),
507 varnames,
508 Geom(),
509 t_new[0], istep, ref_ratio);
510 writeJobInfo(plotfilename);
511#ifdef REMORA_USE_PARTICLES
512 particleData.Checkpoint(plotfilename);
513#endif
514 }
515 }
516 } // end multi-level
517 for (int lev = 0; lev <= finest_level; ++lev) {
518 mask_arrays_for_write(lev, 0.0_rt, (Real) fill_value);
519 }
520}
521
522/**
523 * @param plotfilename name of plotfile to write to
524 * @param nlevels number of levels to write out
525 * @param mf MultiFab of data to write out
526 * @param mf_nd Multifab of nodal data to write out
527 * @param varnames variable names to write out
528 * @param my_geom geometry to use for writing plotfile
529 * @param time time at which to output
530 * @param level_steps vector over level of iterations
531 * @param rr refinement ratio to use for writing plotfile
532 * @param versionName version string for VisIt
533 * @param levelPrefix string to prepend to level number
534 * @param mfPrefix subdirectory for multifab data
535 * @param extra_dirs additional subdirectories within plotfile
536 */
537 void
538 REMORA::WriteMultiLevelPlotfileWithBathymetry (const std::string& plotfilename, int nlevels,
539 const Vector<const MultiFab*>& mf,
540 const Vector<const MultiFab*>& mf_nd,
541 const Vector<const MultiFab*>& mf_u,
542 const Vector<const MultiFab*>& mf_v,
543 const Vector<const MultiFab*>& mf_w,
544 const Vector<std::string>& varnames,
545 const Vector<Geometry>& my_geom,
546 Real time,
547 const Vector<int>& level_steps,
548 const Vector<IntVect>& rr,
549 const std::string &versionName,
550 const std::string &levelPrefix,
551 const std::string &mfPrefix,
552 const Vector<std::string>& extra_dirs) const
553{
554 BL_PROFILE("WriteMultiLevelPlotfileWithBathymetry()");
555
556 BL_ASSERT(nlevels <= mf.size());
557 BL_ASSERT(nlevels <= ref_ratio.size()+1);
558 BL_ASSERT(nlevels <= level_steps.size());
559 BL_ASSERT(mf[0]->nComp() == varnames.size());
560
561 bool callBarrier(false);
562 PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier);
563 if (!extra_dirs.empty()) {
564 for (const auto& d : extra_dirs) {
565 const std::string ed = plotfilename+"/"+d;
566 PreBuildDirectorHierarchy(ed, levelPrefix, nlevels, callBarrier);
567 }
568 }
569 ParallelDescriptor::Barrier();
570
571 if (ParallelDescriptor::MyProc() == ParallelDescriptor::NProcs()-1) {
572 Vector<BoxArray> boxArrays(nlevels);
573 for(int level(0); level < boxArrays.size(); ++level) {
574 boxArrays[level] = mf[level]->boxArray();
575 }
576
577 auto f = [=]() {
578 VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
579 std::string HeaderFileName(plotfilename + "/Header");
580 std::ofstream HeaderFile;
581 HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
582 HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out |
583 std::ofstream::trunc |
584 std::ofstream::binary);
585 if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName);
586 WriteGenericPlotfileHeaderWithBathymetry(HeaderFile, nlevels, boxArrays, varnames,
587 my_geom, time, level_steps, rr, versionName,
588 levelPrefix, mfPrefix);
589 };
590
591 if (AsyncOut::UseAsyncOut()) {
592 AsyncOut::Submit(std::move(f));
593 } else {
594 f();
595 }
596 }
597
598 std::string mf_nodal_prefix = "Nu_nd";
599 std::string mf_uface_prefix = "UFace";
600 std::string mf_vface_prefix = "VFace";
601 std::string mf_wface_prefix = "WFace";
602
603 for (int level = 0; level <= finest_level; ++level)
604 {
605 if (AsyncOut::UseAsyncOut()) {
606 VisMF::AsyncWrite(*mf[level],
607 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix),
608 true);
609 VisMF::AsyncWrite(*mf_nd[level],
610 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix),
611 true);
612 if (plot_staggered_vels) {
613 VisMF::AsyncWrite(*mf_u[level],
614 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_uface_prefix),
615 true);
616 VisMF::AsyncWrite(*mf_v[level],
617 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_vface_prefix),
618 true);
619 VisMF::AsyncWrite(*mf_w[level],
620 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_wface_prefix),
621 true);
622 }
623 } else { const MultiFab* data;
624 std::unique_ptr<MultiFab> mf_tmp;
625 if (mf[level]->nGrowVect() != 0) {
626 mf_tmp = std::make_unique<MultiFab>(mf[level]->boxArray(),
627 mf[level]->DistributionMap(),
628 mf[level]->nComp(), 0, MFInfo(),
629 mf[level]->Factory());
630 MultiFab::Copy(*mf_tmp, *mf[level], 0, 0, mf[level]->nComp(), 0);
631 data = mf_tmp.get();
632 } else {
633 data = mf[level];
634 }
635 VisMF::Write(*data , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix));
636 VisMF::Write(*mf_nd[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix));
637 if (plot_staggered_vels) {
638 VisMF::Write(*mf_u[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_uface_prefix));
639 VisMF::Write(*mf_v[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_vface_prefix));
640 VisMF::Write(*mf_w[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_wface_prefix));
641 }
642 }
643 }
644}
645
646/**
647 * @param HeaderFile output stream for header
648 * @param nlevels number of levels to write out
649 * @param bArray vector over levels of BoxArrays
650 * @param varnames variable names to write out
651 * @param my_geom geometry to use for writing plotfile
652 * @param time time at which to output
653 * @param level_steps vector over level of iterations
654 * @param my_ref_ratio refinement ratio to use for writing plotfile
655 * @param versionName version string for VisIt
656 * @param levelPrefix string to prepend to level number
657 * @param mfPrefix subdirectory for multifab data
658 */
659void
660REMORA::WriteGenericPlotfileHeaderWithBathymetry (std::ostream &HeaderFile,
661 [[maybe_unused]] int nlevels,
662 const Vector<BoxArray> &bArray,
663 const Vector<std::string> &varnames,
664 const Vector<Geometry>& my_geom,
665 Real time,
666 const Vector<int> &level_steps,
667 const Vector<IntVect>& my_ref_ratio,
668 const std::string &versionName,
669 const std::string &levelPrefix,
670 const std::string &mfPrefix) const
671{
672 BL_ASSERT(nlevels <= bArray.size());
673 BL_ASSERT(nlevels <= ref_ratio.size()+1);
674 BL_ASSERT(nlevels <= level_steps.size());
675
676 int num_extra_mfs = 1; // for nodal, which is always on
677 if (plot_staggered_vels) {
678 num_extra_mfs += 3; // for nodal, which is always on
679 }
680
681 HeaderFile.precision(17);
682
683 // ---- this is the generic plot file type name
684 HeaderFile << versionName << '\n';
685
686 HeaderFile << varnames.size() << '\n';
687
688 for (int ivar = 0; ivar < varnames.size(); ++ivar) {
689 HeaderFile << varnames[ivar] << "\n";
690 }
691 HeaderFile << AMREX_SPACEDIM << '\n';
692 HeaderFile << time << '\n';
693 HeaderFile << finest_level << '\n';
694 for (int i = 0; i < AMREX_SPACEDIM; ++i) {
695 HeaderFile << my_geom[0].ProbLo(i) << ' ';
696 }
697 HeaderFile << '\n';
698 for (int i = 0; i < AMREX_SPACEDIM; ++i) {
699 HeaderFile << my_geom[0].ProbHi(i) << ' ';
700 }
701 HeaderFile << '\n';
702 for (int i = 0; i < finest_level; ++i) {
703 HeaderFile << my_ref_ratio[i][0] << ' ';
704 }
705 HeaderFile << '\n';
706 for (int i = 0; i <= finest_level; ++i) {
707 HeaderFile << my_geom[i].Domain() << ' ';
708 }
709 HeaderFile << '\n';
710 for (int i = 0; i <= finest_level; ++i) {
711 HeaderFile << level_steps[i] << ' ';
712 }
713 HeaderFile << '\n';
714 for (int i = 0; i <= finest_level; ++i) {
715 for (int k = 0; k < AMREX_SPACEDIM; ++k) {
716 HeaderFile << my_geom[i].CellSize()[k] << ' ';
717 }
718 HeaderFile << '\n';
719 }
720 HeaderFile << (int) my_geom[0].Coord() << '\n';
721 HeaderFile << "0\n";
722
723 for (int level = 0; level <= finest_level; ++level) {
724 HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n';
725 HeaderFile << level_steps[level] << '\n';
726
727 const IntVect& domain_lo = my_geom[level].Domain().smallEnd();
728 for (int i = 0; i < bArray[level].size(); ++i)
729 {
730 // Need to shift because the RealBox ctor we call takes the
731 // physical location of index (0,0,0). This does not affect
732 // the usual cases where the domain index starts with 0.
733 const Box& b = shift(bArray[level][i], -domain_lo);
734 RealBox loc = RealBox(b, my_geom[level].CellSize(), my_geom[level].ProbLo());
735 for (int n = 0; n < AMREX_SPACEDIM; ++n) {
736 HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n';
737 }
738 }
739
740 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n';
741 }
742 HeaderFile << num_extra_mfs << "\n";
743 HeaderFile << "3" << "\n";
744 HeaderFile << "amrexvec_nu_x" << "\n";
745 HeaderFile << "amrexvec_nu_y" << "\n";
746 HeaderFile << "amrexvec_nu_z" << "\n";
747 std::string mf_nodal_prefix = "Nu_nd";
748 for (int level = 0; level <= finest_level; ++level) {
749 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n';
750 }
751 if (plot_staggered_vels) {
752 HeaderFile << "1" << "\n"; // number of components in the multifab
753 HeaderFile << "u_vel" << "\n";
754 std::string mf_uface_prefix = "UFace";
755 for (int level = 0; level <= finest_level; ++level) {
756 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_uface_prefix) << '\n';
757 }
758 HeaderFile << "1" << "\n";
759 HeaderFile << "v_vel" << "\n";
760 std::string mf_vface_prefix = "VFace";
761 for (int level = 0; level <= finest_level; ++level) {
762 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_vface_prefix) << '\n';
763 }
764 HeaderFile << "1" << "\n";
765 HeaderFile << "w_vel" << "\n";
766 std::string mf_wface_prefix = "WFace";
767 for (int level = 0; level <= finest_level; ++level) {
768 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_wface_prefix) << '\n';
769 }
770 }
771}
772
773/**
774 * @param lev level to mask
775 * @param fill_value fill value to mask with
776 * @param fill_where value at cells where we will apply the mask. This is necessary because rivers
777 */
778void
779REMORA::mask_arrays_for_write(int lev, Real fill_value, Real fill_where) {
780 for (MFIter mfi(*cons_new[lev],false); mfi.isValid(); ++mfi) {
781 Box gbx1 = mfi.growntilebox(IntVect(NGROW+1,NGROW+1,0));
782 Box ubx = mfi.grownnodaltilebox(0,IntVect(NGROW,NGROW,0));
783 Box vbx = mfi.grownnodaltilebox(1,IntVect(NGROW,NGROW,0));
784
785 Array4<Real> const& Zt_avg1 = vec_Zt_avg1[lev]->array(mfi);
786 Array4<Real> const& ubar = vec_ubar[lev]->array(mfi);
787 Array4<Real> const& vbar = vec_vbar[lev]->array(mfi);
788 Array4<Real> const& xvel = xvel_new[lev]->array(mfi);
789 Array4<Real> const& yvel = yvel_new[lev]->array(mfi);
790 Array4<Real> const& temp = cons_new[lev]->array(mfi,Temp_comp);
791 Array4<Real> const& salt = cons_new[lev]->array(mfi,Salt_comp);
792
793 Array4<Real const> const& mskr = vec_mskr[lev]->array(mfi);
794 Array4<Real const> const& msku = vec_msku[lev]->array(mfi);
795 Array4<Real const> const& mskv = vec_mskv[lev]->array(mfi);
796
797 ParallelFor(makeSlab(gbx1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
798 {
799 if (mskr(i,j,0) == 0.0) { // Explicitly compare to 0.0
800 Zt_avg1(i,j,0) = fill_value;
801 }
802 });
803 ParallelFor(gbx1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
804 {
805 if (mskr(i,j,0) == 0.0) { // Explicitly compare to 0.0
806 temp(i,j,k) = fill_value;
807 salt(i,j,k) = fill_value;
808 }
809 });
810 ParallelFor(makeSlab(ubx,2,0), 3, [=] AMREX_GPU_DEVICE (int i, int j, int , int n)
811 {
812 if (msku(i,j,0) == 0.0 && ubar(i,j,0)==fill_where) { // Explicitly compare to 0.0
813 ubar(i,j,0,n) = fill_value;
814 }
815 });
816 ParallelFor(makeSlab(vbx,2,0), 3, [=] AMREX_GPU_DEVICE (int i, int j, int , int n)
817 {
818 if (mskv(i,j,0) == 0.0 && vbar(i,j,0)==fill_where) { // Explicitly compare to 0.0
819 vbar(i,j,0,n) = fill_value;
820 }
821 });
822 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
823 {
824 if (msku(i,j,0) == 0.0 && xvel(i,j,k)==fill_where) { // Explicitly compare to 0.0
825 xvel(i,j,k) = fill_value;
826 }
827 });
828 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
829 {
830 if (mskv(i,j,0) == 0.0 && yvel(i,j,k)==fill_where) { // Explicitly compare to 0.0
831 yvel(i,j,k) = fill_value;
832 }
833 });
834 }
835 Gpu::streamSynchronize();
836}
PlotfileType::amrex
@ amrex
PlotfileType::netcdf
@ netcdf
PlotfileType::hdf5
@ hdf5
Coord
Coord
Coordinates.
Definition REMORA_DataStruct.H:20
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:13
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:8
Salt_comp
#define Salt_comp
Definition REMORA_IndexDefines.H:9
NCONS
#define NCONS
Definition REMORA_IndexDefines.H:11
containerHasElement
bool containerHasElement(const V &iterable, const T &query)
Definition REMORA_Plotfile.cpp:10
null_bc_for_fill
PhysBCFunctNoOp null_bc_for_fill
Definition REMORA_Plotfile.cpp:7
REMORA::plotfile_type
static PlotfileType plotfile_type
Native or NetCDF plotfile output.
Definition REMORA.H:1324
REMORA::cons_names
const amrex::Vector< std::string > cons_names
Names of scalars for plotfile output.
Definition REMORA.H:1274
REMORA::domain_bcs_type
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1179
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:372
REMORA::plot_var_names
amrex::Vector< std::string > plot_var_names
Names of variables to output to AMReX plotfile.
Definition REMORA.H:1272
REMORA::cons_new
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:229
REMORA::vec_mskr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:363
REMORA::writeJobInfo
void writeJobInfo(const std::string &dir) const
Write job info to stdout.
Definition REMORA_writeJobInfo.cpp:7
REMORA::zvel_new
amrex::Vector< amrex::MultiFab * > zvel_new
multilevel data container for current step's z velocities (largely unused; W stored separately)
Definition REMORA.H:235
REMORA::vec_msku
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:365
REMORA::FillPatchNoBC
void FillPatchNoBC(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=true)
Fill a new MultiFab by copying in phi from valid region and filling ghost cells without applying boun...
Definition REMORA_FillPatch.cpp:200
REMORA::setPlotVariables
void setPlotVariables(const std::string &pp_plot_var_names)
Definition REMORA_Plotfile.cpp:18
REMORA::yvel_new
amrex::Vector< amrex::MultiFab * > yvel_new
multilevel data container for current step's y velocities (v in ROMS)
Definition REMORA.H:233
REMORA::xvel_new
amrex::Vector< amrex::MultiFab * > xvel_new
multilevel data container for current step's x velocities (u in ROMS)
Definition REMORA.H:231
REMORA::WriteMultiLevelPlotfileWithBathymetry
void WriteMultiLevelPlotfileWithBathymetry(const std::string &plotfilename, int nlevels, const amrex::Vector< const amrex::MultiFab * > &mf, const amrex::Vector< const amrex::MultiFab * > &mf_nd, const amrex::Vector< const amrex::MultiFab * > &mf_u, const amrex::Vector< const amrex::MultiFab * > &mf_v, const amrex::Vector< const amrex::MultiFab * > &mf_w, const amrex::Vector< std::string > &varnames, const amrex::Vector< amrex::Geometry > &my_geom, amrex::Real time, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &rr, const std::string &versionName="HyperCLaw-V1.1", const std::string &levelPrefix="Level_", const std::string &mfPrefix="Cell", const amrex::Vector< std::string > &extra_dirs=amrex::Vector< std::string >()) const
write out particular data to an AMReX plotfile
Definition REMORA_Plotfile.cpp:538
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:367
REMORA::istep
amrex::Vector< int > istep
which step?
Definition REMORA.H:1153
REMORA::mask_arrays_for_write
void mask_arrays_for_write(int lev, amrex::Real fill_value, amrex::Real fill_where)
Mask data arrays before writing output.
Definition REMORA_Plotfile.cpp:779
REMORA::WriteGenericPlotfileHeaderWithBathymetry
void WriteGenericPlotfileHeaderWithBathymetry(std::ostream &HeaderFile, int nlevels, const amrex::Vector< amrex::BoxArray > &bArray, const amrex::Vector< std::string > &varnames, const amrex::Vector< amrex::Geometry > &my_geom, amrex::Real time, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &rr, const std::string &versionName, const std::string &levelPrefix, const std::string &mfPrefix) const
write out header data for an AMReX plotfile
Definition REMORA_Plotfile.cpp:660
REMORA::file_min_digits
static int file_min_digits
Minimum number of digits in plotfile name or chunked history file.
Definition REMORA.H:1318
REMORA::pp_prefix
std::string pp_prefix
default prefix for input file parameters
Definition REMORA.H:217
REMORA::t_new
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1157
REMORA::plot_staggered_vels
static int plot_staggered_vels
Whether to write the staggered velocities (not averaged to cell centers)
Definition REMORA.H:1321
REMORA::vec_vbar
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:358
REMORA::expand_plotvars_to_unif_rr
bool expand_plotvars_to_unif_rr
whether plotfile variables should be expanded to a uniform refinement ratio
Definition REMORA.H:1283
REMORA::vec_ubar
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:356
REMORA::appendPlotVariables
void appendPlotVariables(const std::string &pp_plot_var_names)
Definition REMORA_Plotfile.cpp:102
REMORA::vec_z_phys_nd
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_phys_nd
z coordinates at psi points (cell nodes)
Definition REMORA.H:293
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:374
REMORA::plot_file_name
std::string plot_file_name
Plotfile prefix.
Definition REMORA.H:1249
REMORA::vec_Zt_avg1
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Zt_avg1
Average of the free surface, zeta (2D)
Definition REMORA.H:296
REMORA::WritePlotFile
void WritePlotFile()
main driver for writing AMReX plotfiles
Definition REMORA_Plotfile.cpp:147
REMORA::derived_names
const amrex::Vector< std::string > derived_names
Names of derived fields for plotfiles.
Definition REMORA.H:1277
derived::remora_dernull
void remora_dernull(const amrex::Box &, amrex::FArrayBox &, int, int, const amrex::FArrayBox &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Geometry &, amrex::Real, const int *, const int)
Definition REMORA_Derive.cpp:9
derived::remora_dervort
void remora_dervort(const amrex::Box &bx, amrex::FArrayBox &derfab, int dcomp, int ncomp, const amrex::FArrayBox &datfab, const amrex::Array4< const amrex::Real > &pm, const amrex::Array4< const amrex::Real > &pn, const amrex::Geometry &, amrex::Real, const int *, const int)
Definition REMORA_Derive.cpp:41