diff --git a/ERF_8H_source.html b/ERF_8H_source.html index 189195121..c1b6ec8e7 100644 --- a/ERF_8H_source.html +++ b/ERF_8H_source.html @@ -1547,7 +1547,7 @@
ERF_WindFarm.H
ERF_WriteBndryPlanes.H
ERF
Definition: ERF.H:124
-
ERF::MakeHorizontalAverages
void MakeHorizontalAverages()
Definition: ERF.cpp:1724
+
ERF::MakeHorizontalAverages
void MakeHorizontalAverages()
Definition: ERF.cpp:1741
ERF::rU_new
amrex::Vector< amrex::MultiFab > rU_new
Definition: ERF.H:758
ERF::last_check_file_step
int last_check_file_step
Definition: ERF.H:910
ERF::walldist
amrex::Vector< std::unique_ptr< amrex::MultiFab > > walldist
Definition: ERF.H:855
@@ -1606,7 +1606,7 @@
ERF::d_havg_qv
amrex::Gpu::DeviceVector< amrex::Real > d_havg_qv
Definition: ERF.H:1158
ERF::column_loc_y
static amrex::Real column_loc_y
Definition: ERF.H:1097
ERF::mg_verbose
static int mg_verbose
Definition: ERF.H:1043
-
ERF::ReadParameters
void ReadParameters()
Definition: ERF.cpp:1377
+
ERF::ReadParameters
void ReadParameters()
Definition: ERF.cpp:1394
ERF::PlotFileVarNames
static amrex::Vector< std::string > PlotFileVarNames(amrex::Vector< std::string > plot_var_names)
Definition: ERF_Plotfile.cpp:174
ERF::InitializeFromFile
void InitializeFromFile()
ERF::physbcs_cons
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_cons > > physbcs_cons
Definition: ERF.H:745
@@ -1642,7 +1642,7 @@
ERF::stop_time
amrex::Real stop_time
Definition: ERF.H:919
ERF::NumSamplePointLogs
AMREX_FORCE_INLINE int NumSamplePointLogs() noexcept
Definition: ERF.H:1259
ERF::h_havg_pressure
amrex::Vector< amrex::Real > h_havg_pressure
Definition: ERF.H:1151
-
ERF::update_diffusive_arrays
void update_diffusive_arrays(int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
Definition: ERF_MakeNewArrays.cpp:382
+
ERF::update_diffusive_arrays
void update_diffusive_arrays(int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
Definition: ERF_MakeNewArrays.cpp:380
ERF::verbose
static int verbose
Definition: ERF.H:1042
ERF::h_havg_qc
amrex::Vector< amrex::Real > h_havg_qc
Definition: ERF.H:1153
ERF::physbcs_w
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_w > > physbcs_w
Definition: ERF.H:748
@@ -1661,7 +1661,7 @@
ERF::Qv_prim
amrex::Vector< std::unique_ptr< amrex::MultiFab > > Qv_prim
Definition: ERF.H:753
ERF::sum_interval
static int sum_interval
Definition: ERF.H:1047
ERF::pert_interval
static int pert_interval
Definition: ERF.H:1048
-
ERF::restart
void restart()
Definition: ERF.cpp:1234
+
ERF::restart
void restart()
Definition: ERF.cpp:1251
ERF::Tau13_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > > Tau13_lev
Definition: ERF.H:810
ERF::SFS_q1fx2_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > > SFS_q1fx2_lev
Definition: ERF.H:822
ERF::ReadCheckpointFileMOST
void ReadCheckpointFileMOST()
Definition: ERF_Checkpoint.cpp:627
@@ -1671,7 +1671,7 @@
ERF::rV_new
amrex::Vector< amrex::MultiFab > rV_new
Definition: ERF.H:760
ERF::PlotFileName
std::string PlotFileName(int lev) const
ERF::write_1D_profiles
void write_1D_profiles(amrex::Real time)
Definition: ERF_Write1DProfiles.cpp:15
-
ERF::initialize_integrator
void initialize_integrator(int lev, amrex::MultiFab &cons_mf, amrex::MultiFab &vel_mf)
Definition: ERF_MakeNewArrays.cpp:590
+
ERF::initialize_integrator
void initialize_integrator(int lev, amrex::MultiFab &cons_mf, amrex::MultiFab &vel_mf)
Definition: ERF_MakeNewArrays.cpp:588
ERF::thin_zforce
amrex::Vector< std::unique_ptr< amrex::MultiFab > > thin_zforce
Definition: ERF.H:905
ERF::init_type
static InitType init_type
Definition: ERF.H:1055
ERF::domain_bcs_type
amrex::Vector< amrex::BCRec > domain_bcs_type
Definition: ERF.H:881
@@ -1702,7 +1702,7 @@
ERF::stretched_dz_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > stretched_dz_d
Definition: ERF.H:863
ERF::t_avg_cnt
amrex::Vector< amrex::Real > t_avg_cnt
Definition: ERF.H:738
ERF::plot_var_names_1
amrex::Vector< std::string > plot_var_names_1
Definition: ERF.H:962
-
ERF::remake_zphys
void remake_zphys(int lev, amrex::Real time, std::unique_ptr< amrex::MultiFab > &temp_zphys_nd)
Definition: ERF_MakeNewArrays.cpp:544
+
ERF::remake_zphys
void remake_zphys(int lev, amrex::Real time, std::unique_ptr< amrex::MultiFab > &temp_zphys_nd)
Definition: ERF_MakeNewArrays.cpp:542
ERF::m_check_int
int m_check_int
Definition: ERF.H:959
ERF::estTimeStep
amrex::Real estTimeStep(int lev, long &dt_fast_ratio) const
Definition: ERF_ComputeTimestep.cpp:55
ERF::sponge_type
static std::string sponge_type
Definition: ERF.H:1059
@@ -1723,11 +1723,11 @@
ERF::output_bndry_planes
static int output_bndry_planes
Definition: ERF.H:1101
ERF::d_v_geos
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > d_v_geos
Definition: ERF.H:1126
ERF::last_plot_file_step_1
int last_plot_file_step_1
Definition: ERF.H:907
-
ERF::update_terrain_arrays
void update_terrain_arrays(int lev)
Definition: ERF_MakeNewArrays.cpp:579
+
ERF::update_terrain_arrays
void update_terrain_arrays(int lev)
Definition: ERF_MakeNewArrays.cpp:577
ERF::nc_bdy_file
static std::string nc_bdy_file
Definition: ERF.H:1069
ERF::metgrid_interp_theta
bool metgrid_interp_theta
Definition: ERF.H:1083
ERF::nc_init_file
static amrex::Vector< amrex::Vector< std::string > > nc_init_file
Definition: ERF.H:1066
-
ERF::init_only
void init_only(int lev, amrex::Real time)
Definition: ERF.cpp:1265
+
ERF::init_only
void init_only(int lev, amrex::Real time)
Definition: ERF.cpp:1282
ERF::input_bndry_planes
static int input_bndry_planes
Definition: ERF.H:1107
ERF::interpolation_type
static StateInterpType interpolation_type
Definition: ERF.H:1056
ERF::d_havg_qc
amrex::Gpu::DeviceVector< amrex::Real > d_havg_qc
Definition: ERF.H:1159
@@ -1757,7 +1757,7 @@
ERF::input_sponge_data
InputSpongeData input_sponge_data
Definition: ERF.H:687
ERF::restart_chkfile
std::string restart_chkfile
Definition: ERF.H:922
ERF::metgrid_use_below_sfc
bool metgrid_use_below_sfc
Definition: ERF.H:1085
-
ERF::init_zphys
void init_zphys(int lev, amrex::Real time)
Definition: ERF_MakeNewArrays.cpp:473
+
ERF::init_zphys
void init_zphys(int lev, amrex::Real time)
Definition: ERF_MakeNewArrays.cpp:471
ERF::sampleptlogname
amrex::Vector< std::string > sampleptlogname
Definition: ERF.H:1378
ERF::xvel_bc_data
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > xvel_bc_data
Definition: ERF.H:690
ERF::InitData_pre
void InitData_pre()
Definition: ERF.cpp:615
@@ -1771,7 +1771,7 @@
ERF::column_loc_x
static amrex::Real column_loc_x
Definition: ERF.H:1096
ERF::ax
amrex::Vector< std::unique_ptr< amrex::MultiFab > > ax
Definition: ERF.H:833
ERF::z_phys_nd
amrex::Vector< std::unique_ptr< amrex::MultiFab > > z_phys_nd
Definition: ERF.H:829
-
ERF::MakeDiagnosticAverage
void MakeDiagnosticAverage(amrex::Vector< amrex::Real > &h_havg, amrex::MultiFab &S, int n)
Definition: ERF.cpp:1830
+
ERF::MakeDiagnosticAverage
void MakeDiagnosticAverage(amrex::Vector< amrex::Real > &h_havg, amrex::MultiFab &S, int n)
Definition: ERF.cpp:1847
ERF::metgrid_proximity
amrex::Real metgrid_proximity
Definition: ERF.H:1088
ERF::h_havg_temperature
amrex::Vector< amrex::Real > h_havg_temperature
Definition: ERF.H:1150
ERF::sampleptlog
amrex::Vector< std::unique_ptr< std::fstream > > sampleptlog
Definition: ERF.H:1377
@@ -1823,7 +1823,7 @@
ERF::ax_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > > ax_src
Definition: ERF.H:839
ERF::input_sponge
void input_sponge(int lev)
Definition: ERF_InitSponge.cpp:17
ERF::d_u_geos
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > d_u_geos
Definition: ERF.H:1123
-
ERF::Define_ERFFillPatchers
void Define_ERFFillPatchers(int lev)
Definition: ERF.cpp:1902
+
ERF::Define_ERFFillPatchers
void Define_ERFFillPatchers(int lev)
Definition: ERF.cpp:1919
ERF::NumDataLogs
AMREX_FORCE_INLINE int NumDataLogs() noexcept
Definition: ERF.H:1245
ERF::turbPert
TurbulentPerturbation turbPert
Definition: ERF.H:1010
ERF::rW_old
amrex::Vector< amrex::MultiFab > rW_old
Definition: ERF.H:761
@@ -1831,7 +1831,7 @@
ERF::FillCoarsePatch
void FillCoarsePatch(int lev, amrex::Real time)
Definition: ERF_FillCoarsePatch.cpp:21
ERF::ClearLevel
void ClearLevel(int lev) override
Definition: ERF_MakeNewLevel.cpp:484
ERF::ref_tags
static amrex::Vector< amrex::AMRErrorTag > ref_tags
Definition: ERF.H:1171
-
ERF::make_physbcs
void make_physbcs(int lev)
Definition: ERF_MakeNewArrays.cpp:612
+
ERF::make_physbcs
void make_physbcs(int lev)
Definition: ERF_MakeNewArrays.cpp:610
ERF::d_w_subsid
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > d_w_subsid
Definition: ERF.H:1120
ERF::FPr_w
amrex::Vector< ERFFillPatcher > FPr_w
Definition: ERF.H:805
ERF::real_set_width
int real_set_width
Definition: ERF.H:1071
@@ -1856,7 +1856,7 @@
ERF::advance_microphysics
void advance_microphysics(int lev, amrex::MultiFab &cons_in, const amrex::Real &dt_advance, const int &iteration, const amrex::Real &time)
Definition: ERF_AdvanceMicrophysics.cpp:5
ERF::Lwave
amrex::Vector< std::unique_ptr< amrex::MultiFab > > Lwave
Definition: ERF.H:870
ERF::project_velocities
void project_velocities(int lev, amrex::Real dt, amrex::Vector< amrex::MultiFab > &vars, amrex::MultiFab &p)
Definition: ERF_PoissonSolve.cpp:10
-
ERF::ParameterSanityChecks
void ParameterSanityChecks()
Definition: ERF.cpp:1658
+
ERF::ParameterSanityChecks
void ParameterSanityChecks()
Definition: ERF.cpp:1675
ERF::derive_stress_profiles
void derive_stress_profiles(amrex::Gpu::HostVector< amrex::Real > &h_avg_tau11, amrex::Gpu::HostVector< amrex::Real > &h_avg_tau12, amrex::Gpu::HostVector< amrex::Real > &h_avg_tau13, amrex::Gpu::HostVector< amrex::Real > &h_avg_tau22, amrex::Gpu::HostVector< amrex::Real > &h_avg_tau23, amrex::Gpu::HostVector< amrex::Real > &h_avg_tau33, amrex::Gpu::HostVector< amrex::Real > &h_avg_hfx3, amrex::Gpu::HostVector< amrex::Real > &h_avg_q1fx3, amrex::Gpu::HostVector< amrex::Real > &h_avg_q2fx3, amrex::Gpu::HostVector< amrex::Real > &h_avg_diss)
Definition: ERF_Write1DProfiles.cpp:482
ERF::Tau12_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > > Tau12_lev
Definition: ERF.H:809
ERF::cf_width
int cf_width
Definition: ERF.H:800
@@ -1864,7 +1864,7 @@
ERF::zflux_imask
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > zflux_imask
Definition: ERF.H:899
ERF::m_expand_plotvars_to_unif_rr
bool m_expand_plotvars_to_unif_rr
Definition: ERF.H:942
ERF::plot_file_on_restart
int plot_file_on_restart
Definition: ERF.H:911
-
ERF::Construct_ERFFillPatchers
void Construct_ERFFillPatchers(int lev)
Definition: ERF.cpp:1876
+
ERF::Construct_ERFFillPatchers
void Construct_ERFFillPatchers(int lev)
Definition: ERF.cpp:1893
ERF::post_update
void post_update(amrex::MultiFab &state_mf, amrex::Real time, const amrex::Geometry &geom)
ERF::num_boxes_at_level
amrex::Vector< int > num_boxes_at_level
Definition: ERF.H:717
ERF::Factory
amrex::FabFactory< amrex::FArrayBox > const & Factory(int lev) const noexcept
Definition: ERF.H:1398
@@ -1907,7 +1907,7 @@
ERF::metgrid_order
int metgrid_order
Definition: ERF.H:1089
ERF::finished_wave
bool finished_wave
Definition: ERF.H:873
ERF::ReadCheckpointFile
void ReadCheckpointFile()
Definition: ERF_Checkpoint.cpp:273
-
ERF::writeNow
bool writeNow(const amrex::Real cur_time, const amrex::Real dt, const int nstep, const int plot_int, const amrex::Real plot_per)
Definition: ERF.cpp:1953
+
ERF::writeNow
bool writeNow(const amrex::Real cur_time, const amrex::Real dt, const int nstep, const int plot_int, const amrex::Real plot_per)
Definition: ERF.cpp:1970
ERF::vars_old
amrex::Vector< amrex::Vector< amrex::MultiFab > > vars_old
Definition: ERF.H:734
ERF::build_fine_mask
amrex::MultiFab & build_fine_mask(int lev)
Definition: ERF_WriteScalarProfiles.cpp:429
ERF::physbcs_base
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_base > > physbcs_base
Definition: ERF.H:749
@@ -1921,7 +1921,7 @@
ERF::Tau22_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > > Tau22_lev
Definition: ERF.H:808
ERF::turbPert_amplitude
void turbPert_amplitude(const int lev)
Definition: ERF_InitTurbPert.cpp:41
ERF::sample_lines
void sample_lines(int lev, amrex::Real time, amrex::IntVect cell, amrex::MultiFab &mf)
Definition: ERF_WriteScalarProfiles.cpp:287
-
ERF::initializeMicrophysics
void initializeMicrophysics(const int &)
Definition: ERF.cpp:1198
+
ERF::initializeMicrophysics
void initializeMicrophysics(const int &)
Definition: ERF.cpp:1215
ERF::plot_lsm
bool plot_lsm
Definition: ERF.H:949
ERF::cons_names
const amrex::Vector< std::string > cons_names
Definition: ERF.H:964
ERF::FillPatch
void FillPatch(int lev, amrex::Real time, const amrex::Vector< amrex::MultiFab * > &mfs_vel, bool cons_only=false)
diff --git a/ERF__NCWpsFile_8H.html b/ERF__NCWpsFile_8H.html index a845446c9..6e8ea8cea 100644 --- a/ERF__NCWpsFile_8H.html +++ b/ERF__NCWpsFile_8H.html @@ -295,55 +295,55 @@

296 {
-
297  int ioproc = amrex::ParallelDescriptor::IOProcessorNumber(); // I/O rank
-
298 
-
299  amrex::Vector<NDArray<float>> nc_arrays(nc_var_names.size());
-
300 
-
301  if (amrex::ParallelDescriptor::IOProcessor())
-
302  {
-
303  ReadNetCDFFile(fname, nc_var_names, nc_arrays);
-
304  }
-
305 
-
306  for (int iv = 0; iv < nc_var_names.size(); iv++)
-
307  {
-
308  FAB tmp;
-
309  if (amrex::ParallelDescriptor::IOProcessor()) {
-
310  fill_fab_from_arrays<FAB,DType>(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp);
-
311  }
-
312 
-
313  int ncomp = tmp.nComp();
-
314  amrex::Box box = tmp.box();
-
315 
-
316  amrex::ParallelDescriptor::Bcast(&box, 1, ioproc);
-
317  amrex::ParallelDescriptor::Bcast(&ncomp, 1, ioproc);
-
318 
-
319  if (!amrex::ParallelDescriptor::IOProcessor()) {
-
320 #ifdef AMREX_USE_GPU
-
321  tmp.resize(box,ncomp,amrex::The_Pinned_Arena());
-
322 #else
-
323  tmp.resize(box,ncomp);
-
324 #endif
-
325  }
+
304 {
+
305  int ioproc = amrex::ParallelDescriptor::IOProcessorNumber(); // I/O rank
+
306 
+
307  amrex::Vector<NDArray<float>> nc_arrays(nc_var_names.size());
+
308 
+
309  if (amrex::ParallelDescriptor::IOProcessor())
+
310  {
+
311  ReadNetCDFFile(fname, nc_var_names, nc_arrays);
+
312  }
+
313 
+
314  for (int iv = 0; iv < nc_var_names.size(); iv++)
+
315  {
+
316  FAB tmp;
+
317  if (amrex::ParallelDescriptor::IOProcessor()) {
+
318  fill_fab_from_arrays<FAB,DType>(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp);
+
319  }
+
320 
+
321  int ncomp = tmp.nComp();
+
322  amrex::Box box = tmp.box();
+
323 
+
324  amrex::ParallelDescriptor::Bcast(&box, 1, ioproc);
+
325  amrex::ParallelDescriptor::Bcast(&ncomp, 1, ioproc);
326 
-
327  amrex::ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc);
-
328 
-
329  // Shift box by the domain lower corner
-
330  amrex::Box fab_bx = tmp.box();
-
331  amrex::Dim3 dom_lb = lbound(domain);
-
332  fab_bx += amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z);
-
333  // fab_vars points to data on device
-
334  fab_vars[iv]->resize(fab_bx,1);
-
335 #ifdef AMREX_USE_GPU
-
336  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
-
337  tmp.dataPtr(), tmp.dataPtr() + tmp.size(),
-
338  fab_vars[iv]->dataPtr());
-
339 #else
-
340  // Provided by BaseFab inheritance through FArrayBox
-
341  fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1);
-
342 #endif
-
343  }
-
344 }
+
327  if (!amrex::ParallelDescriptor::IOProcessor()) {
+
328 #ifdef AMREX_USE_GPU
+
329  tmp.resize(box,ncomp,amrex::The_Pinned_Arena());
+
330 #else
+
331  tmp.resize(box,ncomp);
+
332 #endif
+
333  }
+
334 
+
335  amrex::ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc);
+
336 
+
337  // Shift box by the domain lower corner
+
338  amrex::Box fab_bx = tmp.box();
+
339  amrex::Dim3 dom_lb = lbound(domain);
+
340  fab_bx += amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z);
+
341  // fab_vars points to data on device
+
342  fab_vars[iv]->resize(fab_bx,1);
+
343 #ifdef AMREX_USE_GPU
+
344  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
+
345  tmp.dataPtr(), tmp.dataPtr() + tmp.size(),
+
346  fab_vars[iv]->dataPtr());
+
347 #else
+
348  // Provided by BaseFab inheritance through FArrayBox
+
349  fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1);
+
350 #endif
+
351  }
+
352 }
void ReadNetCDFFile(const std::string &fname, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays)
Definition: ERF_NCWpsFile.H:157
Here is the call graph for this function:
@@ -486,32 +486,40 @@

251  // TODO: The box will only start at (0,0,0) at level 0 -- we need to generalize this
252  amrex::Box my_box(amrex::IntVect(0,0,0), amrex::IntVect(ns3-1,ns2-1,ns1-1));
253 
-
254  if (var_name == "U" || var_name == "UU" ||
-
255  var_name == "MAPFAC_U" || var_name == "MAPFAC_UY") my_box.setType(amrex::IndexType(amrex::IntVect(1,0,0)));
-
256  if (var_name == "V" || var_name == "VV" ||
-
257  var_name == "MAPFAC_V" || var_name == "MAPFAC_VY") my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0)));
-
258  if (var_name == "W" || var_name == "WW") my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
-
259 
-
260  amrex::Arena* Arena_Used = amrex::The_Arena();
-
261 #ifdef AMREX_USE_GPU
-
262  // Make sure temp lives on CPU since nc_arrays lives on CPU only
-
263  Arena_Used = amrex::The_Pinned_Arena();
-
264 #endif
-
265  temp.resize(my_box,1, Arena_Used);
-
266  amrex::Array4<DType> fab_arr = temp.array();
+
254  if (var_name == "PH" || var_name == "PHB") {
+
255  my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
+
256  } else if (var_name == "U" || var_name == "UU" ||
+
257  var_name == "MAPFAC_U" || var_name == "MAPFAC_UY")
+
258  {
+
259  my_box.setType(amrex::IndexType(amrex::IntVect(1,0,0)));
+
260  } else if (var_name == "V" || var_name == "VV" ||
+
261  var_name == "MAPFAC_V" || var_name == "MAPFAC_VY")
+
262  {
+
263  my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0)));
+
264  } else if (var_name == "W" || var_name == "WW") {
+
265  my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
+
266  }
267 
-
268  int ioff = temp.box().smallEnd()[0];
-
269  int joff = temp.box().smallEnd()[1];
-
270 
-
271  auto num_pts = my_box.numPts();
-
272 
-
273  for (int n(0); n < num_pts; ++n) {
-
274  int k = n / (ns2*ns3);
-
275  int j = (n - k*(ns2*ns3)) / ns3 + joff;
-
276  int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff;
-
277  fab_arr(i,j,k,0) = static_cast<DType>(*(nc_arrays[iv].get_data()+n));
-
278  }
-
279 }
+
268  amrex::Arena* Arena_Used = amrex::The_Arena();
+
269 #ifdef AMREX_USE_GPU
+
270  // Make sure temp lives on CPU since nc_arrays lives on CPU only
+
271  Arena_Used = amrex::The_Pinned_Arena();
+
272 #endif
+
273  temp.resize(my_box,1, Arena_Used);
+
274  amrex::Array4<DType> fab_arr = temp.array();
+
275 
+
276  int ioff = temp.box().smallEnd()[0];
+
277  int joff = temp.box().smallEnd()[1];
+
278 
+
279  auto num_pts = my_box.numPts();
+
280 
+
281  for (int n(0); n < num_pts; ++n) {
+
282  int k = n / (ns2*ns3);
+
283  int j = (n - k*(ns2*ns3)) / ns3 + joff;
+
284  int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff;
+
285  fab_arr(i,j,k,0) = static_cast<DType>(*(nc_arrays[iv].get_data()+n));
+
286  }
+
287 }
diff --git a/ERF__NCWpsFile_8H_source.html b/ERF__NCWpsFile_8H_source.html index 6e1b5f520..65f8dbb6c 100644 --- a/ERF__NCWpsFile_8H_source.html +++ b/ERF__NCWpsFile_8H_source.html @@ -339,99 +339,107 @@
251  // TODO: The box will only start at (0,0,0) at level 0 -- we need to generalize this
252  amrex::Box my_box(amrex::IntVect(0,0,0), amrex::IntVect(ns3-1,ns2-1,ns1-1));
253 
-
254  if (var_name == "U" || var_name == "UU" ||
-
255  var_name == "MAPFAC_U" || var_name == "MAPFAC_UY") my_box.setType(amrex::IndexType(amrex::IntVect(1,0,0)));
-
256  if (var_name == "V" || var_name == "VV" ||
-
257  var_name == "MAPFAC_V" || var_name == "MAPFAC_VY") my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0)));
-
258  if (var_name == "W" || var_name == "WW") my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
-
259 
-
260  amrex::Arena* Arena_Used = amrex::The_Arena();
-
261 #ifdef AMREX_USE_GPU
-
262  // Make sure temp lives on CPU since nc_arrays lives on CPU only
-
263  Arena_Used = amrex::The_Pinned_Arena();
-
264 #endif
-
265  temp.resize(my_box,1, Arena_Used);
-
266  amrex::Array4<DType> fab_arr = temp.array();
+
254  if (var_name == "PH" || var_name == "PHB") {
+
255  my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
+
256  } else if (var_name == "U" || var_name == "UU" ||
+
257  var_name == "MAPFAC_U" || var_name == "MAPFAC_UY")
+
258  {
+
259  my_box.setType(amrex::IndexType(amrex::IntVect(1,0,0)));
+
260  } else if (var_name == "V" || var_name == "VV" ||
+
261  var_name == "MAPFAC_V" || var_name == "MAPFAC_VY")
+
262  {
+
263  my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0)));
+
264  } else if (var_name == "W" || var_name == "WW") {
+
265  my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
+
266  }
267 
-
268  int ioff = temp.box().smallEnd()[0];
-
269  int joff = temp.box().smallEnd()[1];
-
270 
-
271  auto num_pts = my_box.numPts();
-
272 
-
273  for (int n(0); n < num_pts; ++n) {
-
274  int k = n / (ns2*ns3);
-
275  int j = (n - k*(ns2*ns3)) / ns3 + joff;
-
276  int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff;
-
277  fab_arr(i,j,k,0) = static_cast<DType>(*(nc_arrays[iv].get_data()+n));
-
278  }
-
279 }
-
280 
-
281 /**
-
282  * Function to read NetCDF variables and fill the corresponding Array4's
-
283  *
-
284  * @param fname Name of the NetCDF file to be read
-
285  * @param nc_var_names Variable names in the NetCDF file
-
286  * @param NC_dim_types NetCDF data dimension types
-
287  * @param fab_vars Fab data we are to fill
-
288  */
-
289 template<class FAB,typename DType>
-
290 void
-
291 BuildFABsFromNetCDFFile (const amrex::Box& domain,
-
292  const std::string &fname,
-
293  amrex::Vector<std::string> nc_var_names,
-
294  amrex::Vector<enum NC_Data_Dims_Type> NC_dim_types,
-
295  amrex::Vector<FAB*> fab_vars)
-
296 {
-
297  int ioproc = amrex::ParallelDescriptor::IOProcessorNumber(); // I/O rank
-
298 
-
299  amrex::Vector<NDArray<float>> nc_arrays(nc_var_names.size());
-
300 
-
301  if (amrex::ParallelDescriptor::IOProcessor())
-
302  {
-
303  ReadNetCDFFile(fname, nc_var_names, nc_arrays);
-
304  }
-
305 
-
306  for (int iv = 0; iv < nc_var_names.size(); iv++)
-
307  {
-
308  FAB tmp;
-
309  if (amrex::ParallelDescriptor::IOProcessor()) {
-
310  fill_fab_from_arrays<FAB,DType>(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp);
-
311  }
-
312 
-
313  int ncomp = tmp.nComp();
-
314  amrex::Box box = tmp.box();
-
315 
-
316  amrex::ParallelDescriptor::Bcast(&box, 1, ioproc);
-
317  amrex::ParallelDescriptor::Bcast(&ncomp, 1, ioproc);
-
318 
-
319  if (!amrex::ParallelDescriptor::IOProcessor()) {
-
320 #ifdef AMREX_USE_GPU
-
321  tmp.resize(box,ncomp,amrex::The_Pinned_Arena());
-
322 #else
-
323  tmp.resize(box,ncomp);
-
324 #endif
-
325  }
+
268  amrex::Arena* Arena_Used = amrex::The_Arena();
+
269 #ifdef AMREX_USE_GPU
+
270  // Make sure temp lives on CPU since nc_arrays lives on CPU only
+
271  Arena_Used = amrex::The_Pinned_Arena();
+
272 #endif
+
273  temp.resize(my_box,1, Arena_Used);
+
274  amrex::Array4<DType> fab_arr = temp.array();
+
275 
+
276  int ioff = temp.box().smallEnd()[0];
+
277  int joff = temp.box().smallEnd()[1];
+
278 
+
279  auto num_pts = my_box.numPts();
+
280 
+
281  for (int n(0); n < num_pts; ++n) {
+
282  int k = n / (ns2*ns3);
+
283  int j = (n - k*(ns2*ns3)) / ns3 + joff;
+
284  int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff;
+
285  fab_arr(i,j,k,0) = static_cast<DType>(*(nc_arrays[iv].get_data()+n));
+
286  }
+
287 }
+
288 
+
289 /**
+
290  * Function to read NetCDF variables and fill the corresponding Array4's
+
291  *
+
292  * @param fname Name of the NetCDF file to be read
+
293  * @param nc_var_names Variable names in the NetCDF file
+
294  * @param NC_dim_types NetCDF data dimension types
+
295  * @param fab_vars Fab data we are to fill
+
296  */
+
297 template<class FAB,typename DType>
+
298 void
+
299 BuildFABsFromNetCDFFile (const amrex::Box& domain,
+
300  const std::string &fname,
+
301  amrex::Vector<std::string> nc_var_names,
+
302  amrex::Vector<enum NC_Data_Dims_Type> NC_dim_types,
+
303  amrex::Vector<FAB*> fab_vars)
+
304 {
+
305  int ioproc = amrex::ParallelDescriptor::IOProcessorNumber(); // I/O rank
+
306 
+
307  amrex::Vector<NDArray<float>> nc_arrays(nc_var_names.size());
+
308 
+
309  if (amrex::ParallelDescriptor::IOProcessor())
+
310  {
+
311  ReadNetCDFFile(fname, nc_var_names, nc_arrays);
+
312  }
+
313 
+
314  for (int iv = 0; iv < nc_var_names.size(); iv++)
+
315  {
+
316  FAB tmp;
+
317  if (amrex::ParallelDescriptor::IOProcessor()) {
+
318  fill_fab_from_arrays<FAB,DType>(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp);
+
319  }
+
320 
+
321  int ncomp = tmp.nComp();
+
322  amrex::Box box = tmp.box();
+
323 
+
324  amrex::ParallelDescriptor::Bcast(&box, 1, ioproc);
+
325  amrex::ParallelDescriptor::Bcast(&ncomp, 1, ioproc);
326 
-
327  amrex::ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc);
-
328 
-
329  // Shift box by the domain lower corner
-
330  amrex::Box fab_bx = tmp.box();
-
331  amrex::Dim3 dom_lb = lbound(domain);
-
332  fab_bx += amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z);
-
333  // fab_vars points to data on device
-
334  fab_vars[iv]->resize(fab_bx,1);
-
335 #ifdef AMREX_USE_GPU
-
336  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
-
337  tmp.dataPtr(), tmp.dataPtr() + tmp.size(),
-
338  fab_vars[iv]->dataPtr());
-
339 #else
-
340  // Provided by BaseFab inheritance through FArrayBox
-
341  fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1);
-
342 #endif
-
343  }
-
344 }
-
345 
-
346 #endif
+
327  if (!amrex::ParallelDescriptor::IOProcessor()) {
+
328 #ifdef AMREX_USE_GPU
+
329  tmp.resize(box,ncomp,amrex::The_Pinned_Arena());
+
330 #else
+
331  tmp.resize(box,ncomp);
+
332 #endif
+
333  }
+
334 
+
335  amrex::ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc);
+
336 
+
337  // Shift box by the domain lower corner
+
338  amrex::Box fab_bx = tmp.box();
+
339  amrex::Dim3 dom_lb = lbound(domain);
+
340  fab_bx += amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z);
+
341  // fab_vars points to data on device
+
342  fab_vars[iv]->resize(fab_bx,1);
+
343 #ifdef AMREX_USE_GPU
+
344  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
+
345  tmp.dataPtr(), tmp.dataPtr() + tmp.size(),
+
346  fab_vars[iv]->dataPtr());
+
347 #else
+
348  // Provided by BaseFab inheritance through FArrayBox
+
349  fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1);
+
350 #endif
+
351  }
+
352 }
+
353 
+
354 #endif
NC_Data_Dims_Type
Definition: ERF_NCWpsFile.H:30
@@ -447,7 +455,7 @@
amrex::Vector< amrex::FArrayBox > PlaneVector
Definition: ERF_NCWpsFile.H:13
void ReadNetCDFFile(const std::string &fname, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays)
Definition: ERF_NCWpsFile.H:157
int BuildFABsFromWRFBdyFile(const std::string &fname, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_xlo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_xhi, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_ylo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_yhi)
-
void BuildFABsFromNetCDFFile(const amrex::Box &domain, const std::string &fname, amrex::Vector< std::string > nc_var_names, amrex::Vector< enum NC_Data_Dims_Type > NC_dim_types, amrex::Vector< FAB * > fab_vars)
Definition: ERF_NCWpsFile.H:291
+
void BuildFABsFromNetCDFFile(const amrex::Box &domain, const std::string &fname, amrex::Vector< std::string > nc_var_names, amrex::Vector< enum NC_Data_Dims_Type > NC_dim_types, amrex::Vector< FAB * > fab_vars)
Definition: ERF_NCWpsFile.H:299
static NCFile open(const std::string &name, const int cmode=NC_NOWRITE)
Definition: ERF_NCInterface.cpp:674
Definition: ERF_NCWpsFile.H:51
bool owned
Definition: ERF_NCWpsFile.H:121
diff --git a/classERF.html b/classERF.html index cb3bb6490..59af4b5a5 100644 --- a/classERF.html +++ b/classERF.html @@ -2919,29 +2919,29 @@

-
1877 {
-
1878  auto& fine_new = vars_new[lev];
-
1879  auto& crse_new = vars_new[lev-1];
-
1880  auto& ba_fine = fine_new[Vars::cons].boxArray();
-
1881  auto& ba_crse = crse_new[Vars::cons].boxArray();
-
1882  auto& dm_fine = fine_new[Vars::cons].DistributionMap();
-
1883  auto& dm_crse = crse_new[Vars::cons].DistributionMap();
-
1884 
-
1885  int ncomp = vars_new[lev][Vars::cons].nComp();
-
1886 
-
1887  FPr_c.emplace_back(ba_fine, dm_fine, geom[lev] ,
-
1888  ba_crse, dm_crse, geom[lev-1],
-
1889  -cf_width, -cf_set_width, ncomp, &cell_cons_interp);
-
1890  FPr_u.emplace_back(convert(ba_fine, IntVect(1,0,0)), dm_fine, geom[lev] ,
-
1891  convert(ba_crse, IntVect(1,0,0)), dm_crse, geom[lev-1],
-
1892  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
-
1893  FPr_v.emplace_back(convert(ba_fine, IntVect(0,1,0)), dm_fine, geom[lev] ,
-
1894  convert(ba_crse, IntVect(0,1,0)), dm_crse, geom[lev-1],
-
1895  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
-
1896  FPr_w.emplace_back(convert(ba_fine, IntVect(0,0,1)), dm_fine, geom[lev] ,
-
1897  convert(ba_crse, IntVect(0,0,1)), dm_crse, geom[lev-1],
-
1898  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
-
1899 }
+
1894 {
+
1895  auto& fine_new = vars_new[lev];
+
1896  auto& crse_new = vars_new[lev-1];
+
1897  auto& ba_fine = fine_new[Vars::cons].boxArray();
+
1898  auto& ba_crse = crse_new[Vars::cons].boxArray();
+
1899  auto& dm_fine = fine_new[Vars::cons].DistributionMap();
+
1900  auto& dm_crse = crse_new[Vars::cons].DistributionMap();
+
1901 
+
1902  int ncomp = vars_new[lev][Vars::cons].nComp();
+
1903 
+
1904  FPr_c.emplace_back(ba_fine, dm_fine, geom[lev] ,
+
1905  ba_crse, dm_crse, geom[lev-1],
+
1906  -cf_width, -cf_set_width, ncomp, &cell_cons_interp);
+
1907  FPr_u.emplace_back(convert(ba_fine, IntVect(1,0,0)), dm_fine, geom[lev] ,
+
1908  convert(ba_crse, IntVect(1,0,0)), dm_crse, geom[lev-1],
+
1909  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
+
1910  FPr_v.emplace_back(convert(ba_fine, IntVect(0,1,0)), dm_fine, geom[lev] ,
+
1911  convert(ba_crse, IntVect(0,1,0)), dm_crse, geom[lev-1],
+
1912  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
+
1913  FPr_w.emplace_back(convert(ba_fine, IntVect(0,0,1)), dm_fine, geom[lev] ,
+
1914  convert(ba_crse, IntVect(0,0,1)), dm_crse, geom[lev-1],
+
1915  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
+
1916 }
int cf_set_width
Definition: ERF.H:801
@@ -3029,29 +3029,29 @@

-
1903 {
-
1904  auto& fine_new = vars_new[lev];
-
1905  auto& crse_new = vars_new[lev-1];
-
1906  auto& ba_fine = fine_new[Vars::cons].boxArray();
-
1907  auto& ba_crse = crse_new[Vars::cons].boxArray();
-
1908  auto& dm_fine = fine_new[Vars::cons].DistributionMap();
-
1909  auto& dm_crse = crse_new[Vars::cons].DistributionMap();
-
1910 
-
1911  int ncomp = fine_new[Vars::cons].nComp();
-
1912 
-
1913  FPr_c[lev-1].Define(ba_fine, dm_fine, geom[lev] ,
-
1914  ba_crse, dm_crse, geom[lev-1],
-
1915  -cf_width, -cf_set_width, ncomp, &cell_cons_interp);
-
1916  FPr_u[lev-1].Define(convert(ba_fine, IntVect(1,0,0)), dm_fine, geom[lev] ,
-
1917  convert(ba_crse, IntVect(1,0,0)), dm_crse, geom[lev-1],
-
1918  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
-
1919  FPr_v[lev-1].Define(convert(ba_fine, IntVect(0,1,0)), dm_fine, geom[lev] ,
-
1920  convert(ba_crse, IntVect(0,1,0)), dm_crse, geom[lev-1],
-
1921  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
-
1922  FPr_w[lev-1].Define(convert(ba_fine, IntVect(0,0,1)), dm_fine, geom[lev] ,
-
1923  convert(ba_crse, IntVect(0,0,1)), dm_crse, geom[lev-1],
-
1924  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
-
1925 }
+
1920 {
+
1921  auto& fine_new = vars_new[lev];
+
1922  auto& crse_new = vars_new[lev-1];
+
1923  auto& ba_fine = fine_new[Vars::cons].boxArray();
+
1924  auto& ba_crse = crse_new[Vars::cons].boxArray();
+
1925  auto& dm_fine = fine_new[Vars::cons].DistributionMap();
+
1926  auto& dm_crse = crse_new[Vars::cons].DistributionMap();
+
1927 
+
1928  int ncomp = fine_new[Vars::cons].nComp();
+
1929 
+
1930  FPr_c[lev-1].Define(ba_fine, dm_fine, geom[lev] ,
+
1931  ba_crse, dm_crse, geom[lev-1],
+
1932  -cf_width, -cf_set_width, ncomp, &cell_cons_interp);
+
1933  FPr_u[lev-1].Define(convert(ba_fine, IntVect(1,0,0)), dm_fine, geom[lev] ,
+
1934  convert(ba_crse, IntVect(1,0,0)), dm_crse, geom[lev-1],
+
1935  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
+
1936  FPr_v[lev-1].Define(convert(ba_fine, IntVect(0,1,0)), dm_fine, geom[lev] ,
+
1937  convert(ba_crse, IntVect(0,1,0)), dm_crse, geom[lev-1],
+
1938  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
+
1939  FPr_w[lev-1].Define(convert(ba_fine, IntVect(0,0,1)), dm_fine, geom[lev] ,
+
1940  convert(ba_crse, IntVect(0,0,1)), dm_crse, geom[lev-1],
+
1941  -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
+
1942 }
@@ -5039,7 +5039,7 @@

amrex::Vector< std::unique_ptr< amrex::MultiFab > > Hwave_onegrid
Definition: ERF.H:871
amrex::Vector< std::unique_ptr< amrex::MultiFab > > thin_yforce
Definition: ERF.H:904
void setPlotVariables(const std::string &pp_plot_var_names, amrex::Vector< std::string > &plot_var_names)
Definition: ERF_Plotfile.cpp:17
-
void ReadParameters()
Definition: ERF.cpp:1377
+
void ReadParameters()
Definition: ERF.cpp:1394
amrex::Vector< std::unique_ptr< amrex::MultiFab > > ay_new
Definition: ERF.H:846
amrex::Vector< std::unique_ptr< amrex::MultiFab > > z_phys_nd_src
Definition: ERF.H:837
amrex::Vector< amrex::MultiFab > base_state_new
Definition: ERF.H:866
@@ -5077,7 +5077,7 @@

amrex::Vector< std::unique_ptr< amrex::MultiFab > > Hwave
Definition: ERF.H:869
amrex::Vector< int > istep
Definition: ERF.H:721
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > xflux_imask
Definition: ERF.H:897
-
void initializeMicrophysics(const int &)
Definition: ERF.cpp:1198
+
void initializeMicrophysics(const int &)
Definition: ERF.cpp:1215
void ReSize(const int &nlev)
Definition: ERF_LandSurface.H:23
const char * buildInfoGetGitHash(int i)
amrex::Real dz0
Definition: ERF_DataStruct.H:637
@@ -5660,7 +5660,7 @@

void WriteCheckpointFile() const
Definition: ERF_Checkpoint.cpp:25
int m_plot_int_2
Definition: ERF.H:944
std::unique_ptr< ReadBndryPlanes > m_r2d
Definition: ERF.H:1164
-
bool writeNow(const amrex::Real cur_time, const amrex::Real dt, const int nstep, const int plot_int, const amrex::Real plot_per)
Definition: ERF.cpp:1953
+
bool writeNow(const amrex::Real cur_time, const amrex::Real dt, const int nstep, const int plot_int, const amrex::Real plot_per)
Definition: ERF.cpp:1970
void timeStep(int lev, amrex::Real time, int iteration)
Definition: ERF_TimeStep.cpp:16
amrex::Real m_plot_per_2
Definition: ERF.H:946
@@ -8145,114 +8145,114 @@

-
1266 {
-
1267  // Map the words in the inputs file to BC types, then translate
-
1268  // those types into what they mean for each variable
-
1269  // This must be called before initHSE (where the base state is initialized)
-
1270  if (lev == 0 && init_type != InitType::Ideal) {
-
1271  init_bcs();
-
1272  }
-
1273 
-
1274  t_new[lev] = time;
-
1275  t_old[lev] = time - 1.e200;
-
1276 
-
1277  auto& lev_new = vars_new[lev];
-
1278  auto& lev_old = vars_old[lev];
-
1279 
-
1280  // Loop over grids at this level to initialize our grid data
-
1281  lev_new[Vars::cons].setVal(0.0); lev_old[Vars::cons].setVal(0.0);
-
1282  lev_new[Vars::xvel].setVal(0.0); lev_old[Vars::xvel].setVal(0.0);
-
1283  lev_new[Vars::yvel].setVal(0.0); lev_old[Vars::yvel].setVal(0.0);
-
1284  lev_new[Vars::zvel].setVal(0.0); lev_old[Vars::zvel].setVal(0.0);
-
1285 
-
1286  // Initialize background flow (optional)
-
1287  if (init_type == InitType::Input_Sounding) {
-
1288  // The base state is initialized by integrating vertically through the
-
1289  // input sounding, if the init_sounding_ideal flag is set; otherwise
-
1290  // it is set by initHSE()
-
1291 
-
1292  // The physbc's need the terrain but are needed for initHSE
-
1293  // We have already made the terrain in the call to init_zphys
-
1294  // in MakeNewLevelFromScratch
-
1295  make_physbcs(lev);
+
1283 {
+
1284  // Map the words in the inputs file to BC types, then translate
+
1285  // those types into what they mean for each variable
+
1286  // This must be called before initHSE (where the base state is initialized)
+
1287  if (lev == 0 && init_type != InitType::Ideal) {
+
1288  init_bcs();
+
1289  }
+
1290 
+
1291  t_new[lev] = time;
+
1292  t_old[lev] = time - 1.e200;
+
1293 
+
1294  auto& lev_new = vars_new[lev];
+
1295  auto& lev_old = vars_old[lev];
1296 
-
1297  // Now init the base state and the data itself
- -
1299 
-
1300  if (init_sounding_ideal) {
-
1301  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(solverChoice.use_gravity,
-
1302  "Gravity should be on to be consistent with sounding initialization.");
-
1303  } else {
-
1304  initHSE();
-
1305  }
-
1306 
-
1307 #ifdef ERF_USE_NETCDF
-
1308  } else if (init_type == InitType::Ideal || init_type == InitType::Real) {
-
1309  // The base state is initialized from WRF wrfinput data, output by
-
1310  // ideal.exe or real.exe
-
1311  init_from_wrfinput(lev);
-
1312 
-
1313  // The physbc's need the terrain but are needed for initHSE
-
1314  if (init_type == InitType::Ideal) {
-
1315  make_physbcs(lev);
-
1316  initHSE(lev);
-
1317  }
-
1318 
-
1319  } else if (init_type == InitType::Metgrid) {
-
1320  // The base state is initialized from data output by WPS metgrid;
-
1321  // we will rebalance after interpolation
-
1322  init_from_metgrid(lev);
-
1323 #endif
-
1324  } else if (init_type == InitType::Uniform) {
-
1325  // Initialize a uniform background field and base state based on the
-
1326  // problem-specified reference density and temperature
-
1327 
-
1328  // The physbc's need the terrain but are needed for initHSE
-
1329  make_physbcs(lev);
-
1330 
-
1331  init_uniform(lev);
-
1332  initHSE(lev);
-
1333  } else {
-
1334  // No background flow initialization specified, initialize the
-
1335  // background field to be equal to the base state, calculated from the
-
1336  // problem-specific erf_init_dens_hse
-
1337 
-
1338  // The bc's need the terrain but are needed for initHSE
-
1339  make_physbcs(lev);
-
1340 
-
1341  // We will initialize the state from the background state so must set that first
-
1342  initHSE(lev);
-
1343  init_from_hse(lev);
-
1344  }
-
1345 
-
1346  // Add problem-specific flow features
-
1347  //
-
1348  // Notes:
-
1349  // - This calls init_custom_pert that is defined for each problem
-
1350  // - This may modify the base state
-
1351  // - The fields set by init_custom_pert are **perturbations** to the
-
1352  // background flow set based on init_type
-
1353  init_custom(lev);
+
1297  // Loop over grids at this level to initialize our grid data
+
1298  lev_new[Vars::cons].setVal(0.0); lev_old[Vars::cons].setVal(0.0);
+
1299  lev_new[Vars::xvel].setVal(0.0); lev_old[Vars::xvel].setVal(0.0);
+
1300  lev_new[Vars::yvel].setVal(0.0); lev_old[Vars::yvel].setVal(0.0);
+
1301  lev_new[Vars::zvel].setVal(0.0); lev_old[Vars::zvel].setVal(0.0);
+
1302 
+
1303  // Initialize background flow (optional)
+
1304  if (init_type == InitType::Input_Sounding) {
+
1305  // The base state is initialized by integrating vertically through the
+
1306  // input sounding, if the init_sounding_ideal flag is set; otherwise
+
1307  // it is set by initHSE()
+
1308 
+
1309  // The physbc's need the terrain but are needed for initHSE
+
1310  // We have already made the terrain in the call to init_zphys
+
1311  // in MakeNewLevelFromScratch
+
1312  make_physbcs(lev);
+
1313 
+
1314  // Now init the base state and the data itself
+ +
1316 
+
1317  if (init_sounding_ideal) {
+
1318  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(solverChoice.use_gravity,
+
1319  "Gravity should be on to be consistent with sounding initialization.");
+
1320  } else {
+
1321  initHSE();
+
1322  }
+
1323 
+
1324 #ifdef ERF_USE_NETCDF
+
1325  } else if (init_type == InitType::Ideal || init_type == InitType::Real) {
+
1326  // The base state is initialized from WRF wrfinput data, output by
+
1327  // ideal.exe or real.exe
+
1328  init_from_wrfinput(lev);
+
1329 
+
1330  // The physbc's need the terrain but are needed for initHSE
+
1331  if (init_type == InitType::Ideal) {
+
1332  make_physbcs(lev);
+
1333  initHSE(lev);
+
1334  }
+
1335 
+
1336  } else if (init_type == InitType::Metgrid) {
+
1337  // The base state is initialized from data output by WPS metgrid;
+
1338  // we will rebalance after interpolation
+
1339  init_from_metgrid(lev);
+
1340 #endif
+
1341  } else if (init_type == InitType::Uniform) {
+
1342  // Initialize a uniform background field and base state based on the
+
1343  // problem-specified reference density and temperature
+
1344 
+
1345  // The physbc's need the terrain but are needed for initHSE
+
1346  make_physbcs(lev);
+
1347 
+
1348  init_uniform(lev);
+
1349  initHSE(lev);
+
1350  } else {
+
1351  // No background flow initialization specified, initialize the
+
1352  // background field to be equal to the base state, calculated from the
+
1353  // problem-specific erf_init_dens_hse
1354 
-
1355  // Ensure that the face-based data are the same on both sides of a periodic domain.
-
1356  // The data associated with the lower grid ID is considered the correct value.
-
1357  lev_new[Vars::xvel].OverrideSync(geom[lev].periodicity());
-
1358  lev_new[Vars::yvel].OverrideSync(geom[lev].periodicity());
-
1359  lev_new[Vars::zvel].OverrideSync(geom[lev].periodicity());
-
1360 
-
1361  if(solverChoice.spongeChoice.sponge_type == "input_sponge"){
-
1362  input_sponge(lev);
-
1363  }
-
1364 
-
1365  // Initialize turbulent perturbation
-
1366  if (solverChoice.pert_type == PerturbationType::Source ||
-
1367  solverChoice.pert_type == PerturbationType::Direct) {
-
1368  if (lev == 0) {
-
1369  turbPert_update(lev, 0.);
-
1370  turbPert_amplitude(lev);
-
1371  }
-
1372  }
-
1373 }
+
1355  // The bc's need the terrain but are needed for initHSE
+
1356  make_physbcs(lev);
+
1357 
+
1358  // We will initialize the state from the background state so must set that first
+
1359  initHSE(lev);
+
1360  init_from_hse(lev);
+
1361  }
+
1362 
+
1363  // Add problem-specific flow features
+
1364  //
+
1365  // Notes:
+
1366  // - This calls init_custom_pert that is defined for each problem
+
1367  // - This may modify the base state
+
1368  // - The fields set by init_custom_pert are **perturbations** to the
+
1369  // background flow set based on init_type
+
1370  init_custom(lev);
+
1371 
+
1372  // Ensure that the face-based data are the same on both sides of a periodic domain.
+
1373  // The data associated with the lower grid ID is considered the correct value.
+
1374  lev_new[Vars::xvel].OverrideSync(geom[lev].periodicity());
+
1375  lev_new[Vars::yvel].OverrideSync(geom[lev].periodicity());
+
1376  lev_new[Vars::zvel].OverrideSync(geom[lev].periodicity());
+
1377 
+
1378  if(solverChoice.spongeChoice.sponge_type == "input_sponge"){
+
1379  input_sponge(lev);
+
1380  }
+
1381 
+
1382  // Initialize turbulent perturbation
+
1383  if (solverChoice.pert_type == PerturbationType::Source ||
+
1384  solverChoice.pert_type == PerturbationType::Direct) {
+
1385  if (lev == 0) {
+
1386  turbPert_update(lev, 0.);
+
1387  turbPert_amplitude(lev);
+
1388  }
+
1389  }
+
1390 }
void init_from_input_sounding(int lev)
Definition: ERF_InitFromInputSounding.cpp:50
static InitType init_type
Definition: ERF.H:1055
void init_custom(int lev)
Definition: ERF_InitCustom.cpp:26
@@ -8261,7 +8261,7 @@

void initHSE()
Initialize HSE.
Definition: ERF_Init1D.cpp:130

void turbPert_update(const int lev, const amrex::Real dt)
Definition: ERF_InitTurbPert.cpp:12
void input_sponge(int lev)
Definition: ERF_InitSponge.cpp:17
-
void make_physbcs(int lev)
Definition: ERF_MakeNewArrays.cpp:612
+
void make_physbcs(int lev)
Definition: ERF_MakeNewArrays.cpp:610
void init_uniform(int lev)
Definition: ERF_InitUniform.cpp:17
void turbPert_amplitude(const int lev)
Definition: ERF_InitTurbPert.cpp:41
bool use_gravity
Definition: ERF_DataStruct.H:614
@@ -8353,337 +8353,335 @@

46  SolverChoice::mesh_type == MeshType::VariableDz) {

47  z_phys_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
48 
-
49  amrex::Print() << "WE ARE LIVE WITH ZPHYS! " << std::endl;
-
50 
-
51  if (solverChoice.terrain_type == TerrainType::MovingFittedMesh)
-
52  {
-
53  detJ_cc_new[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
-
54  detJ_cc_src[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
-
55 
-
56  ax_src[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
-
57  ay_src[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
-
58  az_src[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
-
59 
-
60  ax_new[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
-
61  ay_new[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
-
62  az_new[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
-
63 
-
64  z_t_rk[lev] = std::make_unique<MultiFab>( convert(ba, IntVect(0,0,1)), dm, 1, 1 );
-
65  }
-
66 
-
67  BoxArray ba_nd(ba);
-
68  ba_nd.surroundingNodes();
-
69 
-
70  // We need this to be one greater than the ghost cells to handle levels > 0
- -
72  tmp_zphys_nd = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
-
73 
-
74  if (solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
-
75  z_phys_nd_new[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
-
76  z_phys_nd_src[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
-
77  }
-
78 
-
79  } else {
-
80  z_phys_nd[lev] = nullptr;
-
81  z_phys_cc[lev] = nullptr;
-
82 
-
83  z_phys_nd_new[lev] = nullptr;
-
84  detJ_cc_new[lev] = nullptr;
-
85 
-
86  z_phys_nd_src[lev] = nullptr;
-
87  detJ_cc_src[lev] = nullptr;
-
88 
-
89  z_t_rk[lev] = nullptr;
-
90  }
-
91 
-
92  if (SolverChoice::terrain_type == TerrainType::ImmersedForcing)
-
93  {
-
94  terrain_blanking[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
-
95  terrain_blanking[lev]->setVal(1.0);
-
96  }
-
97 
-
98  // We use these area arrays regardless of terrain, EB or none of the above
-
99  detJ_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
-
100  ax[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
-
101  ay[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
-
102  az[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
-
103 
-
104  detJ_cc[lev]->setVal(1.0);
-
105  ax[lev]->setVal(1.0);
-
106  ay[lev]->setVal(1.0);
-
107  az[lev]->setVal(1.0);
-
108 
+
49  if (solverChoice.terrain_type == TerrainType::MovingFittedMesh)
+
50  {
+
51  detJ_cc_new[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
+
52  detJ_cc_src[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
+
53 
+
54  ax_src[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
+
55  ay_src[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
+
56  az_src[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
+
57 
+
58  ax_new[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
+
59  ay_new[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
+
60  az_new[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
+
61 
+
62  z_t_rk[lev] = std::make_unique<MultiFab>( convert(ba, IntVect(0,0,1)), dm, 1, 1 );
+
63  }
+
64 
+
65  BoxArray ba_nd(ba);
+
66  ba_nd.surroundingNodes();
+
67 
+
68  // We need this to be one greater than the ghost cells to handle levels > 0
+ +
70  tmp_zphys_nd = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
+
71 
+
72  if (solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
+
73  z_phys_nd_new[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
+
74  z_phys_nd_src[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
+
75  }
+
76 
+
77  } else {
+
78  z_phys_nd[lev] = nullptr;
+
79  z_phys_cc[lev] = nullptr;
+
80 
+
81  z_phys_nd_new[lev] = nullptr;
+
82  detJ_cc_new[lev] = nullptr;
+
83 
+
84  z_phys_nd_src[lev] = nullptr;
+
85  detJ_cc_src[lev] = nullptr;
+
86 
+
87  z_t_rk[lev] = nullptr;
+
88  }
+
89 
+
90  if (SolverChoice::terrain_type == TerrainType::ImmersedForcing)
+
91  {
+
92  terrain_blanking[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
+
93  terrain_blanking[lev]->setVal(1.0);
+
94  }
+
95 
+
96  // We use these area arrays regardless of terrain, EB or none of the above
+
97  detJ_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
+
98  ax[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
+
99  ay[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
+
100  az[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
+
101 
+
102  detJ_cc[lev]->setVal(1.0);
+
103  ax[lev]->setVal(1.0);
+
104  ay[lev]->setVal(1.0);
+
105  az[lev]->setVal(1.0);
+
106 
+
107  // ********************************************************************************************
+
108  // Create wall distance array for RANS modeling
109  // ********************************************************************************************
-
110  // Create wall distance array for RANS modeling
-
111  // ********************************************************************************************
-
112  if (solverChoice.turbChoice[lev].rans_type != RANSType::None) {
-
113  walldist[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
-
114  walldist[lev]->setVal(1e23);
-
115  } else {
-
116  walldist[lev] = nullptr;
-
117  }
-
118 
+
110  if (solverChoice.turbChoice[lev].rans_type != RANSType::None) {
+
111  walldist[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
+
112  walldist[lev]->setVal(1e23);
+
113  } else {
+
114  walldist[lev] = nullptr;
+
115  }
+
116 
+
117  // ********************************************************************************************
+
118  // These are the persistent containers for the old and new data
119  // ********************************************************************************************
-
120  // These are the persistent containers for the old and new data
-
121  // ********************************************************************************************
-
122  int ncomp;
-
123  if (lev > 0) {
-
124  ncomp = vars_new[lev-1][Vars::cons].nComp();
-
125  } else {
-
126  int n_qstate = micro->Get_Qstate_Size();
-
127  ncomp = NDRY + NSCALARS + n_qstate;
-
128  }
-
129 
-
130  // ********************************************************************************************
-
131  // The number of ghost cells for density must be 1 greater than that for velocity
-
132  // so that we can go back in forth between velocity and momentum on all faces
-
133  // ********************************************************************************************
- - -
136 
+
120  int ncomp;
+
121  if (lev > 0) {
+
122  ncomp = vars_new[lev-1][Vars::cons].nComp();
+
123  } else {
+
124  int n_qstate = micro->Get_Qstate_Size();
+
125  ncomp = NDRY + NSCALARS + n_qstate;
+
126  }
+
127 
+
128  // ********************************************************************************************
+
129  // The number of ghost cells for density must be 1 greater than that for velocity
+
130  // so that we can go back in forth between velocity and momentum on all faces
+
131  // ********************************************************************************************
+ + +
134 
+
135  // ********************************************************************************************
+
136  // New solution data containers
137  // ********************************************************************************************
-
138  // New solution data containers
-
139  // ********************************************************************************************
-
140  lev_new[Vars::cons].define(ba, dm, ncomp, ngrow_state);
-
141  lev_old[Vars::cons].define(ba, dm, ncomp, ngrow_state);
-
142 
-
143  lev_new[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
-
144  lev_old[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
-
145 
-
146  lev_new[Vars::yvel].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
-
147  lev_old[Vars::yvel].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
-
148 
-
149  // Note that we need the ghost cells in the z-direction if we are doing any
-
150  // kind of domain decomposition in the vertical (at level 0 or above)
-
151  lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
-
152  lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
-
153 
-
154  if (solverChoice.anelastic[lev] == 1) {
-
155  pp_inc[lev].define(ba, dm, 1, 1);
-
156  pp_inc[lev].setVal(0.0);
-
157  }
-
158 
+
138  lev_new[Vars::cons].define(ba, dm, ncomp, ngrow_state);
+
139  lev_old[Vars::cons].define(ba, dm, ncomp, ngrow_state);
+
140 
+
141  lev_new[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
+
142  lev_old[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
+
143 
+
144  lev_new[Vars::yvel].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
+
145  lev_old[Vars::yvel].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
+
146 
+
147  // Note that we need the ghost cells in the z-direction if we are doing any
+
148  // kind of domain decomposition in the vertical (at level 0 or above)
+
149  lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
+
150  lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
+
151 
+
152  if (solverChoice.anelastic[lev] == 1) {
+
153  pp_inc[lev].define(ba, dm, 1, 1);
+
154  pp_inc[lev].setVal(0.0);
+
155  }
+
156 
+
157  // ********************************************************************************************
+
158  // These are just used for scratch in the time integrator but we might as well define them here
159  // ********************************************************************************************
-
160  // These are just used for scratch in the time integrator but we might as well define them here
-
161  // ********************************************************************************************
-
162  rU_old[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
-
163  rU_new[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
-
164 
-
165  rV_old[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
-
166  rV_new[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
-
167 
-
168  rW_old[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
-
169  rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
-
170 
-
171  if (lev > 0) {
-
172  //xmom_crse_rhs[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, IntVect{0});
-
173  //ymom_crse_rhs[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, IntVect{0});
-
174  zmom_crse_rhs[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect{0});
-
175  }
-
176 
-
177  // We do this here just so they won't be undefined in the initial FillPatch
-
178  rU_old[lev].setVal(1.2e21);
-
179  rV_old[lev].setVal(3.4e22);
-
180  rW_old[lev].setVal(5.6e23);
-
181  rU_new[lev].setVal(1.2e21);
-
182  rV_new[lev].setVal(3.4e22);
-
183  rW_new[lev].setVal(5.6e23);
-
184 
+
160  rU_old[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
+
161  rU_new[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
+
162 
+
163  rV_old[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
+
164  rV_new[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
+
165 
+
166  rW_old[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
+
167  rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
+
168 
+
169  if (lev > 0) {
+
170  //xmom_crse_rhs[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, IntVect{0});
+
171  //ymom_crse_rhs[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, IntVect{0});
+
172  zmom_crse_rhs[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect{0});
+
173  }
+
174 
+
175  // We do this here just so they won't be undefined in the initial FillPatch
+
176  rU_old[lev].setVal(1.2e21);
+
177  rV_old[lev].setVal(3.4e22);
+
178  rW_old[lev].setVal(5.6e23);
+
179  rU_new[lev].setVal(1.2e21);
+
180  rV_new[lev].setVal(3.4e22);
+
181  rW_new[lev].setVal(5.6e23);
+
182 
+
183  // ********************************************************************************************
+
184  // These are just time averaged fields for diagnostics
185  // ********************************************************************************************
-
186  // These are just time averaged fields for diagnostics
-
187  // ********************************************************************************************
-
188 
-
189  // NOTE: We are not completing a fillpach call on the time averaged data;
-
190  // which would copy on intersection and interpolate from coarse.
-
191  // Therefore, we are restarting the averaging when the ba changes,
-
192  // this may give poor statistics for dynamic mesh refinement.
-
193  vel_t_avg[lev] = nullptr;
- -
195  vel_t_avg[lev] = std::make_unique<MultiFab>(ba, dm, 4, 0); // Each vel comp and the mag
-
196  vel_t_avg[lev]->setVal(0.0);
-
197  t_avg_cnt[lev] = 0.0;
-
198  }
-
199 
+
186 
+
187  // NOTE: We are not completing a fillpach call on the time averaged data;
+
188  // which would copy on intersection and interpolate from coarse.
+
189  // Therefore, we are restarting the averaging when the ba changes,
+
190  // this may give poor statistics for dynamic mesh refinement.
+
191  vel_t_avg[lev] = nullptr;
+ +
193  vel_t_avg[lev] = std::make_unique<MultiFab>(ba, dm, 4, 0); // Each vel comp and the mag
+
194  vel_t_avg[lev]->setVal(0.0);
+
195  t_avg_cnt[lev] = 0.0;
+
196  }
+
197 
+
198  // ********************************************************************************************
+
199  // Initialize flux registers whenever we create/re-create a level
200  // ********************************************************************************************
-
201  // Initialize flux registers whenever we create/re-create a level
-
202  // ********************************************************************************************
-
203  if (solverChoice.coupling_type == CouplingType::TwoWay) {
-
204  if (lev == 0) {
-
205  advflux_reg[0] = nullptr;
-
206  } else {
-
207  int ncomp_reflux = vars_new[0][Vars::cons].nComp();
-
208  advflux_reg[lev] = new YAFluxRegister(ba , grids[lev-1],
-
209  dm , dmap[lev-1],
-
210  geom[lev], geom[lev-1],
-
211  ref_ratio[lev-1], lev, ncomp_reflux);
-
212  }
-
213  }
-
214 
+
201  if (solverChoice.coupling_type == CouplingType::TwoWay) {
+
202  if (lev == 0) {
+
203  advflux_reg[0] = nullptr;
+
204  } else {
+
205  int ncomp_reflux = vars_new[0][Vars::cons].nComp();
+
206  advflux_reg[lev] = new YAFluxRegister(ba , grids[lev-1],
+
207  dm , dmap[lev-1],
+
208  geom[lev], geom[lev-1],
+
209  ref_ratio[lev-1], lev, ncomp_reflux);
+
210  }
+
211  }
+
212 
+
213  // ********************************************************************************************
+
214  // Define Theta_prim storage if using MOST BC
215  // ********************************************************************************************
-
216  // Define Theta_prim storage if using MOST BC
-
217  // ********************************************************************************************
-
218  if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) {
-
219  Theta_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
-
220  if (solverChoice.moisture_type != MoistureType::None) {
-
221  Qv_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
-
222  Qr_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
-
223  } else {
-
224  Qv_prim[lev] = nullptr;
-
225  Qr_prim[lev] = nullptr;
-
226  }
-
227  } else {
-
228  Theta_prim[lev] = nullptr;
-
229  Qv_prim[lev] = nullptr;
-
230  Qr_prim[lev] = nullptr;
-
231  }
-
232 
+
216  if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) {
+
217  Theta_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
+
218  if (solverChoice.moisture_type != MoistureType::None) {
+
219  Qv_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
+
220  Qr_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
+
221  } else {
+
222  Qv_prim[lev] = nullptr;
+
223  Qr_prim[lev] = nullptr;
+
224  }
+
225  } else {
+
226  Theta_prim[lev] = nullptr;
+
227  Qv_prim[lev] = nullptr;
+
228  Qr_prim[lev] = nullptr;
+
229  }
+
230 
+
231  // ********************************************************************************************
+
232  // Map factors
233  // ********************************************************************************************
-
234  // Map factors
-
235  // ********************************************************************************************
-
236  BoxList bl2d_mf = ba.boxList();
-
237  for (auto& b : bl2d_mf) {
-
238  b.setRange(2,0);
-
239  }
-
240  BoxArray ba2d_mf(std::move(bl2d_mf));
-
241 
-
242  mapfac_m[lev] = std::make_unique<MultiFab>(ba2d_mf,dm,1,3);
-
243  mapfac_u[lev] = std::make_unique<MultiFab>(convert(ba2d_mf,IntVect(1,0,0)),dm,1,3);
-
244  mapfac_v[lev] = std::make_unique<MultiFab>(convert(ba2d_mf,IntVect(0,1,0)),dm,1,3);
- -
246  mapfac_m[lev]->setVal(0.5);
-
247  mapfac_u[lev]->setVal(0.5);
-
248  mapfac_v[lev]->setVal(0.5);
-
249  }
-
250  else {
-
251  mapfac_m[lev]->setVal(1.);
-
252  mapfac_u[lev]->setVal(1.);
-
253  mapfac_v[lev]->setVal(1.);
-
254  }
-
255 
-
256 #if defined(ERF_USE_WINDFARM)
+
234  BoxList bl2d_mf = ba.boxList();
+
235  for (auto& b : bl2d_mf) {
+
236  b.setRange(2,0);
+
237  }
+
238  BoxArray ba2d_mf(std::move(bl2d_mf));
+
239 
+
240  mapfac_m[lev] = std::make_unique<MultiFab>(ba2d_mf,dm,1,3);
+
241  mapfac_u[lev] = std::make_unique<MultiFab>(convert(ba2d_mf,IntVect(1,0,0)),dm,1,3);
+
242  mapfac_v[lev] = std::make_unique<MultiFab>(convert(ba2d_mf,IntVect(0,1,0)),dm,1,3);
+ +
244  mapfac_m[lev]->setVal(0.5);
+
245  mapfac_u[lev]->setVal(0.5);
+
246  mapfac_v[lev]->setVal(0.5);
+
247  }
+
248  else {
+
249  mapfac_m[lev]->setVal(1.);
+
250  mapfac_u[lev]->setVal(1.);
+
251  mapfac_v[lev]->setVal(1.);
+
252  }
+
253 
+
254 #if defined(ERF_USE_WINDFARM)
+
255  //*********************************************************
+
256  // Variables for Fitch model for windfarm parametrization
257  //*********************************************************
-
258  // Variables for Fitch model for windfarm parametrization
-
259  //*********************************************************
-
260  if (solverChoice.windfarm_type == WindFarmType::Fitch){
-
261  vars_windfarm[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt
-
262  }
-
263  if (solverChoice.windfarm_type == WindFarmType::EWP){
-
264  vars_windfarm[lev].define(ba, dm, 3, ngrow_state); // dudt, dvdt, dTKEdt
-
265  }
-
266  if (solverChoice.windfarm_type == WindFarmType::SimpleAD) {
-
267  vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt
-
268  }
-
269  if (solverChoice.windfarm_type == WindFarmType::GeneralAD) {
-
270  vars_windfarm[lev].define(ba, dm, 3, ngrow_state);// dudt, dvdt, dwdt
-
271  }
-
272  Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell
-
273  SMark[lev].define(ba, dm, 2, 1); // Free stream velocity/source term
-
274  // sampling marker in a cell - 2 components
-
275 #endif
-
276 
-
277 
-
278 #ifdef ERF_USE_WW3_COUPLING
-
279  // create a new BoxArray and DistributionMapping for a MultiFab with 1 box
-
280  BoxArray ba_onegrid(geom[lev].Domain());
-
281  BoxList bl2d_onegrid = ba_onegrid.boxList();
-
282  for (auto& b : bl2d_onegrid) {
-
283  b.setRange(2,0);
-
284  }
-
285  BoxArray ba2d_onegrid(std::move(bl2d_onegrid));
-
286  Vector<int> pmap;
-
287  pmap.resize(1);
-
288  pmap[0]=0;
-
289  DistributionMapping dm_onegrid(ba2d_onegrid);
-
290  dm_onegrid.define(pmap);
-
291 
-
292  Hwave_onegrid[lev] = std::make_unique<MultiFab>(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0));
-
293  Lwave_onegrid[lev] = std::make_unique<MultiFab>(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0));
-
294 
-
295  BoxList bl2d_wave = ba.boxList();
-
296  for (auto& b : bl2d_wave) {
-
297  b.setRange(2,0);
-
298  }
-
299  BoxArray ba2d_wave(std::move(bl2d_wave));
-
300 
-
301  Hwave[lev] = std::make_unique<MultiFab>(ba2d_wave,dm,1,IntVect(3,3,0));
-
302  Lwave[lev] = std::make_unique<MultiFab>(ba2d_wave,dm,1,IntVect(3,3,0));
-
303 
-
304  std::cout<<ba_onegrid<<std::endl;
-
305  std::cout<<ba2d_onegrid<<std::endl;
-
306  std::cout<<dm_onegrid<<std::endl;
-
307 #endif
-
308 
-
309 
-
310 #if defined(ERF_USE_RRTMGP)
+
258  if (solverChoice.windfarm_type == WindFarmType::Fitch){
+
259  vars_windfarm[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt
+
260  }
+
261  if (solverChoice.windfarm_type == WindFarmType::EWP){
+
262  vars_windfarm[lev].define(ba, dm, 3, ngrow_state); // dudt, dvdt, dTKEdt
+
263  }
+
264  if (solverChoice.windfarm_type == WindFarmType::SimpleAD) {
+
265  vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt
+
266  }
+
267  if (solverChoice.windfarm_type == WindFarmType::GeneralAD) {
+
268  vars_windfarm[lev].define(ba, dm, 3, ngrow_state);// dudt, dvdt, dwdt
+
269  }
+
270  Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell
+
271  SMark[lev].define(ba, dm, 2, 1); // Free stream velocity/source term
+
272  // sampling marker in a cell - 2 components
+
273 #endif
+
274 
+
275 
+
276 #ifdef ERF_USE_WW3_COUPLING
+
277  // create a new BoxArray and DistributionMapping for a MultiFab with 1 box
+
278  BoxArray ba_onegrid(geom[lev].Domain());
+
279  BoxList bl2d_onegrid = ba_onegrid.boxList();
+
280  for (auto& b : bl2d_onegrid) {
+
281  b.setRange(2,0);
+
282  }
+
283  BoxArray ba2d_onegrid(std::move(bl2d_onegrid));
+
284  Vector<int> pmap;
+
285  pmap.resize(1);
+
286  pmap[0]=0;
+
287  DistributionMapping dm_onegrid(ba2d_onegrid);
+
288  dm_onegrid.define(pmap);
+
289 
+
290  Hwave_onegrid[lev] = std::make_unique<MultiFab>(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0));
+
291  Lwave_onegrid[lev] = std::make_unique<MultiFab>(ba2d_onegrid,dm_onegrid,1,IntVect(1,1,0));
+
292 
+
293  BoxList bl2d_wave = ba.boxList();
+
294  for (auto& b : bl2d_wave) {
+
295  b.setRange(2,0);
+
296  }
+
297  BoxArray ba2d_wave(std::move(bl2d_wave));
+
298 
+
299  Hwave[lev] = std::make_unique<MultiFab>(ba2d_wave,dm,1,IntVect(3,3,0));
+
300  Lwave[lev] = std::make_unique<MultiFab>(ba2d_wave,dm,1,IntVect(3,3,0));
+
301 
+
302  std::cout<<ba_onegrid<<std::endl;
+
303  std::cout<<ba2d_onegrid<<std::endl;
+
304  std::cout<<dm_onegrid<<std::endl;
+
305 #endif
+
306 
+
307 
+
308 #if defined(ERF_USE_RRTMGP)
+
309  //*********************************************************
+
310  // Radiation heating source terms
311  //*********************************************************
-
312  // Radiation heating source terms
-
313  //*********************************************************
-
314  qheating_rates[lev] = std::make_unique<MultiFab>(ba, dm, 2, ngrow_state);
-
315  qheating_rates[lev]->setVal(0.);
-
316 
+
312  qheating_rates[lev] = std::make_unique<MultiFab>(ba, dm, 2, ngrow_state);
+
313  qheating_rates[lev]->setVal(0.);
+
314 
+
315  //*********************************************************
+
316  // Radiation fluxes for coupling to LSM
317  //*********************************************************
-
318  // Radiation fluxes for coupling to LSM
-
319  //*********************************************************
-
320 
-
321  // NOTE: Finer levels do not need to coincide with the bottom domain boundary
-
322  // at k=0. We make slabs here with the kmin for a given box. Therefore,
-
323  // care must be taken before applying these fluxes to an LSM model. For
-
324 
-
325  // Radiative fluxes for LSM
-
326  if (solverChoice.lsm_type != LandSurfaceType::None)
-
327  {
-
328  BoxList m_bl = ba.boxList();
-
329  for (auto& b : m_bl) {
-
330  int kmin = b.smallEnd(2);
-
331  b.setRange(2,kmin);
-
332  }
-
333  BoxArray m_ba(std::move(m_bl));
-
334 
-
335  sw_lw_fluxes[lev] = std::make_unique<MultiFab>(m_ba, dm, 5, ngrow_state); // SW direct (2), SW diffuse (2), LW
-
336  solar_zenith[lev] = std::make_unique<MultiFab>(m_ba, dm, 2, ngrow_state);
-
337 
-
338  sw_lw_fluxes[lev]->setVal(0.);
-
339  solar_zenith[lev]->setVal(0.);
-
340  }
-
341 #endif
-
342 
+
318 
+
319  // NOTE: Finer levels do not need to coincide with the bottom domain boundary
+
320  // at k=0. We make slabs here with the kmin for a given box. Therefore,
+
321  // care must be taken before applying these fluxes to an LSM model. For
+
322 
+
323  // Radiative fluxes for LSM
+
324  if (solverChoice.lsm_type != LandSurfaceType::None)
+
325  {
+
326  BoxList m_bl = ba.boxList();
+
327  for (auto& b : m_bl) {
+
328  int kmin = b.smallEnd(2);
+
329  b.setRange(2,kmin);
+
330  }
+
331  BoxArray m_ba(std::move(m_bl));
+
332 
+
333  sw_lw_fluxes[lev] = std::make_unique<MultiFab>(m_ba, dm, 5, ngrow_state); // SW direct (2), SW diffuse (2), LW
+
334  solar_zenith[lev] = std::make_unique<MultiFab>(m_ba, dm, 2, ngrow_state);
+
335 
+
336  sw_lw_fluxes[lev]->setVal(0.);
+
337  solar_zenith[lev]->setVal(0.);
+
338  }
+
339 #endif
+
340 
+
341  //*********************************************************
+
342  // Turbulent perturbation region initialization
343  //*********************************************************
-
344  // Turbulent perturbation region initialization
-
345  //*********************************************************
-
346  // TODO: Test perturbation on multiple levels
-
347  if (solverChoice.pert_type == PerturbationType::Source ||
-
348  solverChoice.pert_type == PerturbationType::Direct)
-
349  {
-
350  if (lev == 0) {
-
351  turbPert.init_tpi(lev, geom[lev].Domain().bigEnd(), geom[lev].CellSizeArray(), ba, dm, ngrow_state);
-
352  }
-
353  }
-
354 
-
355  //
-
356  // Define the land mask here and set it to all land by default
-
357  // NOTE: the logic below will BREAK if we have any grids not touching the bottom boundary
-
358  //
-
359  {
-
360  lmask_lev[lev].resize(1);
-
361  auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0;
-
362  BoxList bl2d_mask = ba.boxList();
-
363  for (auto& b : bl2d_mask) {
-
364  b.setRange(2,0);
-
365  }
-
366  BoxArray ba2d_mask(std::move(bl2d_mask));
-
367  lmask_lev[lev][0] = std::make_unique<iMultiFab>(ba2d_mask,dm,1,ngv);
-
368  lmask_lev[lev][0]->setVal(1);
-
369  lmask_lev[lev][0]->FillBoundary(geom[lev].periodicity());
-
370  }
-
371 
-
372  // Read in tables needed for windfarm simulations
-
373  // fill in Nturb multifab - number of turbines in each mesh cell
-
374  // write out the vtk files for wind turbine location and/or
-
375  // actuator disks
-
376  #ifdef ERF_USE_WINDFARM
-
377  //init_windfarm(lev);
-
378  #endif
-
379 }
+
344  // TODO: Test perturbation on multiple levels
+
345  if (solverChoice.pert_type == PerturbationType::Source ||
+
346  solverChoice.pert_type == PerturbationType::Direct)
+
347  {
+
348  if (lev == 0) {
+
349  turbPert.init_tpi(lev, geom[lev].Domain().bigEnd(), geom[lev].CellSizeArray(), ba, dm, ngrow_state);
+
350  }
+
351  }
+
352 
+
353  //
+
354  // Define the land mask here and set it to all land by default
+
355  // NOTE: the logic below will BREAK if we have any grids not touching the bottom boundary
+
356  //
+
357  {
+
358  lmask_lev[lev].resize(1);
+
359  auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0;
+
360  BoxList bl2d_mask = ba.boxList();
+
361  for (auto& b : bl2d_mask) {
+
362  b.setRange(2,0);
+
363  }
+
364  BoxArray ba2d_mask(std::move(bl2d_mask));
+
365  lmask_lev[lev][0] = std::make_unique<iMultiFab>(ba2d_mask,dm,1,ngv);
+
366  lmask_lev[lev][0]->setVal(1);
+
367  lmask_lev[lev][0]->FillBoundary(geom[lev].periodicity());
+
368  }
+
369 
+
370  // Read in tables needed for windfarm simulations
+
371  // fill in Nturb multifab - number of turbines in each mesh cell
+
372  // write out the vtk files for wind turbine location and/or
+
373  // actuator disks
+
374  #ifdef ERF_USE_WINDFARM
+
375  //init_windfarm(lev);
+
376  #endif
+
377 }
#define NDRY
Definition: ERF_IndexDefines.H:13
static AMREX_FORCE_INLINE int ComputeGhostCells(const AdvChoice &advChoice, bool use_num_diff)
Definition: ERF.H:1183
@ num_comps
Definition: ERF_IndexDefines.H:67
@@ -8906,74 +8904,74 @@

-
474 {
-
475  if (init_type != InitType::Real && init_type != InitType::Metgrid)
-
476  {
-
477  if (lev > 0 && z_phys_nd[lev]) {
-
478  //
-
479  // First interpolate from coarser level if there is one
-
480  // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab
-
481  // have been pre-filled - this includes ghost cells both inside and outside
-
482  // the domain
-
483  //
-
484  InterpFromCoarseLevel(*z_phys_nd[lev], z_phys_nd[lev]->nGrowVect(),
-
485  IntVect(0,0,0), // do not fill ghost cells outside the domain
-
486  *z_phys_nd[lev-1], 0, 0, 1,
-
487  geom[lev-1], geom[lev],
-
488  refRatio(lev-1), &node_bilinear_interp,
- -
490  }
-
491 
- -
493  Box bx(surroundingNodes(Geom(0).Domain())); bx.grow(ngrow);
-
494  FArrayBox terrain_fab(makeSlab(bx,2,0),1);
-
495 
-
496  //
-
497  // If we are using fitted mesh then we use the surface as defined above
-
498  // If we are not using fitted mesh but are using z_levels, we still need z_phys (for now)
-
499  // but we need to use a flat terrain for the mesh itself (the EB data has already been made
-
500  // from the correct terrain)
-
501  //
-
502  if (solverChoice.terrain_type != TerrainType::StaticFittedMesh &&
-
503  solverChoice.terrain_type != TerrainType::MovingFittedMesh) {
-
504  terrain_fab.template setVal<RunOn::Device>(0.0);
-
505  } else {
+
472 {
+
473  if (init_type != InitType::Real && init_type != InitType::Metgrid)
+
474  {
+
475  if (lev > 0 && z_phys_nd[lev]) {
+
476  //
+
477  // First interpolate from coarser level if there is one
+
478  // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab
+
479  // have been pre-filled - this includes ghost cells both inside and outside
+
480  // the domain
+
481  //
+
482  InterpFromCoarseLevel(*z_phys_nd[lev], z_phys_nd[lev]->nGrowVect(),
+
483  IntVect(0,0,0), // do not fill ghost cells outside the domain
+
484  *z_phys_nd[lev-1], 0, 0, 1,
+
485  geom[lev-1], geom[lev],
+
486  refRatio(lev-1), &node_bilinear_interp,
+ +
488  }
+
489 
+ +
491  Box bx(surroundingNodes(Geom(0).Domain())); bx.grow(ngrow);
+
492  FArrayBox terrain_fab(makeSlab(bx,2,0),1);
+
493 
+
494  //
+
495  // If we are using fitted mesh then we use the surface as defined above
+
496  // If we are not using fitted mesh but are using z_levels, we still need z_phys (for now)
+
497  // but we need to use a flat terrain for the mesh itself (the EB data has already been made
+
498  // from the correct terrain)
+
499  //
+
500  if (solverChoice.terrain_type != TerrainType::StaticFittedMesh &&
+
501  solverChoice.terrain_type != TerrainType::MovingFittedMesh) {
+
502  terrain_fab.template setVal<RunOn::Device>(0.0);
+
503  } else {
+
504  //
+
505  // Fill the values of the terrain height at k=0 only
506  //
-
507  // Fill the values of the terrain height at k=0 only
-
508  //
-
509  prob->init_terrain_surface(geom[lev],terrain_fab,time);
-
510  }
-
511 
-
512  if (z_phys_nd[lev]) { // Has this been allocated?
-
513  for (MFIter mfi(*z_phys_nd[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
-
514  {
-
515  Box isect = terrain_fab.box() & (*z_phys_nd[lev])[mfi].box();
-
516  (*z_phys_nd[lev])[mfi].template copy<RunOn::Device>(terrain_fab,isect,0,isect,0,1);
-
517  }
- -
519  }
-
520 
-
521  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
-
522  terrain_blanking[lev]->setVal(1.0);
-
523  MultiFab::Subtract(*terrain_blanking[lev], EBFactory(lev).getVolFrac(), 0, 0, 1, 0);
-
524  }
-
525 
-
526  if (lev == 0 && z_phys_nd[0]) {
-
527  Real zmax = z_phys_nd[0]->max(0,0,false);
-
528  Real rel_diff = (zmax - zlevels_stag[0][zlevels_stag[0].size()-1]) / zmax;
-
529  if (rel_diff < 1.e-8) {
-
530  amrex::Print() << "max of zphys_nd " << zmax << std::endl;
-
531  amrex::Print() << "max of zlevels " << zlevels_stag[0][zlevels_stag[0].size()-1] << std::endl;
-
532  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rel_diff < 1.e-8, "Terrain is taller than domain top!");
-
533  }
-
534  } // lev == 0
-
535 
-
536  if (z_phys_nd[lev]) {
-
537  z_phys_nd[lev]->FillBoundary(geom[lev].periodicity());
-
538  }
-
539 
-
540  } // init_type
-
541 }
+
507  prob->init_terrain_surface(geom[lev],terrain_fab,time);
+
508  }
+
509 
+
510  if (z_phys_nd[lev]) { // Has this been allocated?
+
511  for (MFIter mfi(*z_phys_nd[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
+
512  {
+
513  Box isect = terrain_fab.box() & (*z_phys_nd[lev])[mfi].box();
+
514  (*z_phys_nd[lev])[mfi].template copy<RunOn::Device>(terrain_fab,isect,0,isect,0,1);
+
515  }
+ +
517  }
+
518 
+
519  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+
520  terrain_blanking[lev]->setVal(1.0);
+
521  MultiFab::Subtract(*terrain_blanking[lev], EBFactory(lev).getVolFrac(), 0, 0, 1, 0);
+
522  }
+
523 
+
524  if (lev == 0 && z_phys_nd[0]) {
+
525  Real zmax = z_phys_nd[0]->max(0,0,false);
+
526  Real rel_diff = (zmax - zlevels_stag[0][zlevels_stag[0].size()-1]) / zmax;
+
527  if (rel_diff < 1.e-8) {
+
528  amrex::Print() << "max of zphys_nd " << zmax << std::endl;
+
529  amrex::Print() << "max of zlevels " << zlevels_stag[0][zlevels_stag[0].size()-1] << std::endl;
+
530  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rel_diff < 1.e-8, "Terrain is taller than domain top!");
+
531  }
+
532  } // lev == 0
+
533 
+
534  if (z_phys_nd[lev]) {
+
535  z_phys_nd[lev]->FillBoundary(geom[lev].periodicity());
+
536  }
+
537 
+
538  } // init_type
+
539 }
void make_terrain_fitted_coords(int lev, const Geometry &geom, MultiFab &z_phys_nd, Vector< Real > const &z_levels_h, GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &phys_bc_type)
Definition: ERF_TerrainMetrics.cpp:118
Here is the call graph for this function:
@@ -9255,326 +9253,343 @@

872  auto& lev_new = vars_new[lev];

873  auto& lev_old = vars_old[lev];
874 
-
875  int ncomp = lev_new[Vars::cons].nComp();
-
876 
+
875  // ***************************************************************************
+
876  // Physical bc's at domain boundary
877  // ***************************************************************************
-
878  // Physical bc's at domain boundary
-
879  // ***************************************************************************
-
880  IntVect ngvect_cons = vars_new[lev][Vars::cons].nGrowVect();
-
881  IntVect ngvect_vels = vars_new[lev][Vars::xvel].nGrowVect();
-
882 
-
883  (*physbcs_cons[lev])(lev_new[Vars::cons],0,ncomp,ngvect_cons,t_new[lev],BCVars::cons_bc,true);
-
884  ( *physbcs_u[lev])(lev_new[Vars::xvel],0,1 ,ngvect_vels,t_new[lev],BCVars::xvel_bc,true);
-
885  ( *physbcs_v[lev])(lev_new[Vars::yvel],0,1 ,ngvect_vels,t_new[lev],BCVars::yvel_bc,true);
-
886  ( *physbcs_w[lev])(lev_new[Vars::zvel],lev_new[Vars::xvel],lev_new[Vars::yvel],
-
887  ngvect_vels,t_new[lev],BCVars::zvel_bc,true);
-
888 
-
889  MultiFab::Copy(lev_old[Vars::cons],lev_new[Vars::cons],0,0,ncomp,lev_new[Vars::cons].nGrowVect());
-
890  MultiFab::Copy(lev_old[Vars::xvel],lev_new[Vars::xvel],0,0, 1,lev_new[Vars::xvel].nGrowVect());
-
891  MultiFab::Copy(lev_old[Vars::yvel],lev_new[Vars::yvel],0,0, 1,lev_new[Vars::yvel].nGrowVect());
-
892  MultiFab::Copy(lev_old[Vars::zvel],lev_new[Vars::zvel],0,0, 1,lev_new[Vars::zvel].nGrowVect());
-
893  }
-
894 
-
895  // Compute the minimum dz in the domain at each level (to be used for setting the timestep)
-
896  dz_min.resize(max_level+1);
-
897  for (int lev = 0; lev <= finest_level; ++lev)
-
898  {
-
899  dz_min[lev] = geom[lev].CellSize(2);
-
900  if ( SolverChoice::mesh_type != MeshType::ConstantDz ) {
-
901  dz_min[lev] *= (*detJ_cc[lev]).min(0);
-
902  }
-
903  }
-
904 
-
905  ComputeDt();
-
906 
-
907  // Fill ghost cells/faces
-
908  for (int lev = 0; lev <= finest_level; ++lev)
-
909  {
-
910  if (lev > 0 && cf_width >= 0) {
- -
912  }
-
913 
-
914  auto& lev_new = vars_new[lev];
-
915 
-
916  //
-
917  // Fill boundary conditions -- not sure why we need this here
-
918  //
-
919  bool fillset = false;
-
920  if (lev == 0) {
-
921  FillPatch(lev, t_new[lev],
-
922  {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]});
-
923  } else {
-
924  FillPatch(lev, t_new[lev],
-
925  {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]},
-
926  {&lev_new[Vars::cons],&rU_new[lev],&rV_new[lev],&rW_new[lev]},
-
927  base_state[lev], base_state[lev],
-
928  fillset);
+
878  IntVect ngvect_cons = vars_new[lev][Vars::cons].nGrowVect();
+
879  IntVect ngvect_vels = vars_new[lev][Vars::xvel].nGrowVect();
+
880 
+
881  int ncomp_cons = lev_new[Vars::cons].nComp();
+
882  bool do_fb = true;
+
883 
+
884 #ifdef ERF_USE_NETCDF
+
885  // We call this here because it is an ERF routine
+
886  if (use_real_bcs && (lev==0)) {
+
887  int icomp_cons = 0;
+
888  bool cons_only = false;
+
889  Vector<MultiFab*> mfs_vec = {&lev_new[Vars::cons],&lev_new[Vars::xvel],
+
890  &lev_new[Vars::yvel],&lev_new[Vars::zvel]};
+
891  fill_from_realbdy(mfs_vec,t_new[lev],cons_only,icomp_cons,
+
892  ncomp_cons,ngvect_cons,ngvect_vels);
+
893  do_fb = false;
+
894  }
+
895 #endif
+
896 
+
897  (*physbcs_cons[lev])(lev_new[Vars::cons],0,ncomp_cons,
+
898  ngvect_cons,t_new[lev],BCVars::cons_bc,do_fb);
+
899  ( *physbcs_u[lev])(lev_new[Vars::xvel],0,1 ,
+
900  ngvect_vels,t_new[lev],BCVars::xvel_bc,do_fb);
+
901  ( *physbcs_v[lev])(lev_new[Vars::yvel],0,1 ,
+
902  ngvect_vels,t_new[lev],BCVars::yvel_bc,do_fb);
+
903  ( *physbcs_w[lev])(lev_new[Vars::zvel],lev_new[Vars::xvel],lev_new[Vars::yvel],
+
904  ngvect_vels,t_new[lev],BCVars::zvel_bc,do_fb);
+
905 
+
906  MultiFab::Copy(lev_old[Vars::cons],lev_new[Vars::cons],0,0,ncomp_cons,lev_new[Vars::cons].nGrowVect());
+
907  MultiFab::Copy(lev_old[Vars::xvel],lev_new[Vars::xvel],0,0, 1,lev_new[Vars::xvel].nGrowVect());
+
908  MultiFab::Copy(lev_old[Vars::yvel],lev_new[Vars::yvel],0,0, 1,lev_new[Vars::yvel].nGrowVect());
+
909  MultiFab::Copy(lev_old[Vars::zvel],lev_new[Vars::zvel],0,0, 1,lev_new[Vars::zvel].nGrowVect());
+
910  }
+
911 
+
912  // Compute the minimum dz in the domain at each level (to be used for setting the timestep)
+
913  dz_min.resize(max_level+1);
+
914  for (int lev = 0; lev <= finest_level; ++lev)
+
915  {
+
916  dz_min[lev] = geom[lev].CellSize(2);
+
917  if ( SolverChoice::mesh_type != MeshType::ConstantDz ) {
+
918  dz_min[lev] *= (*detJ_cc[lev]).min(0);
+
919  }
+
920  }
+
921 
+
922  ComputeDt();
+
923 
+
924  // Fill ghost cells/faces
+
925  for (int lev = 0; lev <= finest_level; ++lev)
+
926  {
+
927  if (lev > 0 && cf_width >= 0) {
+
929  }
930 
-
931  //
-
932  // We do this here to make sure level (lev-1) boundary conditions are filled
-
933  // before we interpolate to level (lev) ghost cells
-
934  //
-
935  if (lev < finest_level) {
-
936  auto& lev_old = vars_old[lev];
-
937  MultiFab::Copy(lev_old[Vars::cons],lev_new[Vars::cons],0,0,lev_old[Vars::cons].nComp(),lev_old[Vars::cons].nGrowVect());
-
938  MultiFab::Copy(lev_old[Vars::xvel],lev_new[Vars::xvel],0,0,lev_old[Vars::xvel].nComp(),lev_old[Vars::xvel].nGrowVect());
-
939  MultiFab::Copy(lev_old[Vars::yvel],lev_new[Vars::yvel],0,0,lev_old[Vars::yvel].nComp(),lev_old[Vars::yvel].nGrowVect());
-
940  MultiFab::Copy(lev_old[Vars::zvel],lev_new[Vars::zvel],0,0,lev_old[Vars::zvel].nComp(),lev_old[Vars::zvel].nGrowVect());
-
941  }
-
942 
-
943  //
-
944  // We fill the ghost cell values of the base state in case it wasn't done in the initialization
-
945  //
-
946  base_state[lev].FillBoundary(geom[lev].periodicity());
+
931  auto& lev_new = vars_new[lev];
+
932 
+
933  //
+
934  // Fill boundary conditions -- not sure why we need this here
+
935  //
+
936  bool fillset = false;
+
937  if (lev == 0) {
+
938  FillPatch(lev, t_new[lev],
+
939  {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]});
+
940  } else {
+
941  FillPatch(lev, t_new[lev],
+
942  {&lev_new[Vars::cons],&lev_new[Vars::xvel],&lev_new[Vars::yvel],&lev_new[Vars::zvel]},
+
943  {&lev_new[Vars::cons],&rU_new[lev],&rV_new[lev],&rW_new[lev]},
+
944  base_state[lev], base_state[lev],
+
945  fillset);
+
946  }
947 
-
948  // For moving terrain only
-
949  if (solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
-
950  MultiFab::Copy(base_state_new[lev],base_state[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
-
951  base_state_new[lev].FillBoundary(geom[lev].periodicity());
-
952  }
-
953 
-
954  }
-
955 
-
956  // Allow idealized cases over water, used to set lmask
-
957  ParmParse pp("erf");
-
958  int is_land;
-
959  for (int lev = 0; lev <= finest_level; ++lev)
-
960  {
-
961  if (pp.query("is_land", is_land, lev)) {
-
962  if (is_land == 1) {
-
963  amrex::Print() << "Level " << lev << " is land" << std::endl;
-
964  } else if (is_land == 0) {
-
965  amrex::Print() << "Level " << lev << " is water" << std::endl;
-
966  } else {
-
967  Error("is_land should be 0 or 1");
-
968  }
-
969  lmask_lev[lev][0]->setVal(is_land);
-
970  lmask_lev[lev][0]->FillBoundary(geom[lev].periodicity());
-
971  }
-
972  }
-
973 
-
974 #ifdef ERF_USE_WW3_COUPLING
-
975  int lev = 0;
-
976  amrex::Print() << " About to call send_to_ww3 from ERF.cpp" << std::endl;
-
977  send_to_ww3(lev);
-
978  amrex::Print() << " About to call read_waves from ERF.cpp" << std::endl;
-
979  read_waves(lev);
-
980  // send_to_ww3(lev);
-
981 #endif
-
982 
-
983  // Configure ABLMost params if used MostWall boundary condition
-
984  // NOTE: we must set up the MOST routine after calling FillPatch
-
985  // in order to have lateral ghost cells filled (MOST + terrain interp).
-
986  // FillPatch does not call MOST, FillIntermediatePatch does.
-
987  if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST)
-
988  {
-
989  bool use_exp_most = solverChoice.use_explicit_most;
-
990  bool use_rot_most = solverChoice.use_rotate_most;
-
991  if (use_exp_most) {
-
992  Print() << "Using MOST with explicitly included surface stresses" << std::endl;
-
993  if (use_rot_most) {
-
994  Print() << "Using MOST with surface stress rotations" << std::endl;
-
995  }
-
996  }
-
997 
-
998  m_most = std::make_unique<ABLMost>(geom, use_exp_most, use_rot_most,
- - - -
1002 #ifdef ERF_USE_NETCDF
-
1003  ,start_bdy_time, bdy_time_interval
-
1004 #endif
-
1005  );
-
1006 
-
1007 
-
1008  if (restart_chkfile != "") {
-
1009  // Update surface fields if needed
- -
1011  }
-
1012 
-
1013  // We now configure ABLMost params here so that we can print the averages at t=0
-
1014  // Note we don't fill ghost cells here because this is just for diagnostics
-
1015  for (int lev = 0; lev <= finest_level; ++lev)
-
1016  {
-
1017  Real time = t_new[lev];
-
1018  IntVect ng = Theta_prim[lev]->nGrowVect();
-
1019 
-
1020  MultiFab::Copy( *Theta_prim[lev], vars_new[lev][Vars::cons], RhoTheta_comp, 0, 1, ng);
-
1021  MultiFab::Divide(*Theta_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);
-
1022 
-
1023  if (solverChoice.moisture_type != MoistureType::None) {
-
1024  ng = Qv_prim[lev]->nGrowVect();
-
1025 
-
1026  MultiFab::Copy( *Qv_prim[lev], vars_new[lev][Vars::cons], RhoQ1_comp, 0, 1, ng);
-
1027  MultiFab::Divide(*Qv_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);
-
1028 
-
1029  int rhoqr_comp = solverChoice.RhoQr_comp;
-
1030  if (rhoqr_comp > -1) {
-
1031  MultiFab::Copy( *Qr_prim[lev], vars_new[lev][Vars::cons], rhoqr_comp, 0, 1, ng);
-
1032  MultiFab::Divide(*Qr_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);
-
1033  } else {
-
1034  Qr_prim[lev]->setVal(0.0);
-
1035  }
-
1036  }
-
1037  m_most->update_mac_ptrs(lev, vars_new, Theta_prim, Qv_prim, Qr_prim);
-
1038 
-
1039  if (restart_chkfile == "") {
-
1040  // Only do this if starting from scratch; if restarting, then
-
1041  // we don't want to call update_fluxes multiple times because
-
1042  // it will change u* and theta* from their previous values
-
1043  m_most->update_pblh(lev, vars_new, z_phys_cc[lev].get(),
- - - -
1047  m_most->update_fluxes(lev, time);
-
1048  }
-
1049  }
-
1050  }
-
1051 
-
1052  // Update micro vars before first plot file
-
1053  if (solverChoice.moisture_type != MoistureType::None) {
-
1054  for (int lev = 0; lev <= finest_level; ++lev) micro->Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
-
1055  }
-
1056 
-
1057  // Fill time averaged velocities before first plot file
-
1058  if (solverChoice.time_avg_vel) {
-
1059  for (int lev = 0; lev <= finest_level; ++lev) {
-
1060  Time_Avg_Vel_atCC(dt[lev], t_avg_cnt[lev], vel_t_avg[lev].get(),
-
1061  vars_new[lev][Vars::xvel],
-
1062  vars_new[lev][Vars::yvel],
-
1063  vars_new[lev][Vars::zvel]);
-
1064  }
-
1065  }
-
1066 
-
1067  // check for additional plotting variables that are available after particle containers
-
1068  // are setup.
-
1069  const std::string& pv1 = "plot_vars_1"; appendPlotVariables(pv1,plot_var_names_1);
-
1070  const std::string& pv2 = "plot_vars_2"; appendPlotVariables(pv2,plot_var_names_2);
-
1071 
-
1072  if ( restart_chkfile.empty() && (m_check_int > 0 || m_check_per > 0.) )
-
1073  {
-
1074 #ifdef ERF_USE_NETCDF
-
1075  if (check_type == "netcdf") {
-
1076  WriteNCCheckpointFile();
-
1077  }
-
1078 #endif
-
1079  if (check_type == "native") {
- +
948  //
+
949  // We do this here to make sure level (lev-1) boundary conditions are filled
+
950  // before we interpolate to level (lev) ghost cells
+
951  //
+
952  if (lev < finest_level) {
+
953  auto& lev_old = vars_old[lev];
+
954  MultiFab::Copy(lev_old[Vars::cons],lev_new[Vars::cons],0,0,lev_old[Vars::cons].nComp(),lev_old[Vars::cons].nGrowVect());
+
955  MultiFab::Copy(lev_old[Vars::xvel],lev_new[Vars::xvel],0,0,lev_old[Vars::xvel].nComp(),lev_old[Vars::xvel].nGrowVect());
+
956  MultiFab::Copy(lev_old[Vars::yvel],lev_new[Vars::yvel],0,0,lev_old[Vars::yvel].nComp(),lev_old[Vars::yvel].nGrowVect());
+
957  MultiFab::Copy(lev_old[Vars::zvel],lev_new[Vars::zvel],0,0,lev_old[Vars::zvel].nComp(),lev_old[Vars::zvel].nGrowVect());
+
958  }
+
959 
+
960  //
+
961  // We fill the ghost cell values of the base state in case it wasn't done in the initialization
+
962  //
+
963  base_state[lev].FillBoundary(geom[lev].periodicity());
+
964 
+
965  // For moving terrain only
+
966  if (solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
+
967  MultiFab::Copy(base_state_new[lev],base_state[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
+
968  base_state_new[lev].FillBoundary(geom[lev].periodicity());
+
969  }
+
970 
+
971  }
+
972 
+
973  // Allow idealized cases over water, used to set lmask
+
974  ParmParse pp("erf");
+
975  int is_land;
+
976  for (int lev = 0; lev <= finest_level; ++lev)
+
977  {
+
978  if (pp.query("is_land", is_land, lev)) {
+
979  if (is_land == 1) {
+
980  amrex::Print() << "Level " << lev << " is land" << std::endl;
+
981  } else if (is_land == 0) {
+
982  amrex::Print() << "Level " << lev << " is water" << std::endl;
+
983  } else {
+
984  Error("is_land should be 0 or 1");
+
985  }
+
986  lmask_lev[lev][0]->setVal(is_land);
+
987  lmask_lev[lev][0]->FillBoundary(geom[lev].periodicity());
+
988  }
+
989  }
+
990 
+
991 #ifdef ERF_USE_WW3_COUPLING
+
992  int lev = 0;
+
993  amrex::Print() << " About to call send_to_ww3 from ERF.cpp" << std::endl;
+
994  send_to_ww3(lev);
+
995  amrex::Print() << " About to call read_waves from ERF.cpp" << std::endl;
+
996  read_waves(lev);
+
997  // send_to_ww3(lev);
+
998 #endif
+
999 
+
1000  // Configure ABLMost params if used MostWall boundary condition
+
1001  // NOTE: we must set up the MOST routine after calling FillPatch
+
1002  // in order to have lateral ghost cells filled (MOST + terrain interp).
+
1003  // FillPatch does not call MOST, FillIntermediatePatch does.
+
1004  if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST)
+
1005  {
+
1006  bool use_exp_most = solverChoice.use_explicit_most;
+
1007  bool use_rot_most = solverChoice.use_rotate_most;
+
1008  if (use_exp_most) {
+
1009  Print() << "Using MOST with explicitly included surface stresses" << std::endl;
+
1010  if (use_rot_most) {
+
1011  Print() << "Using MOST with surface stress rotations" << std::endl;
+
1012  }
+
1013  }
+
1014 
+
1015  m_most = std::make_unique<ABLMost>(geom, use_exp_most, use_rot_most,
+ + + +
1019 #ifdef ERF_USE_NETCDF
+
1020  ,start_bdy_time, bdy_time_interval
+
1021 #endif
+
1022  );
+
1023 
+
1024 
+
1025  if (restart_chkfile != "") {
+
1026  // Update surface fields if needed
+ +
1028  }
+
1029 
+
1030  // We now configure ABLMost params here so that we can print the averages at t=0
+
1031  // Note we don't fill ghost cells here because this is just for diagnostics
+
1032  for (int lev = 0; lev <= finest_level; ++lev)
+
1033  {
+
1034  Real time = t_new[lev];
+
1035  IntVect ng = Theta_prim[lev]->nGrowVect();
+
1036 
+
1037  MultiFab::Copy( *Theta_prim[lev], vars_new[lev][Vars::cons], RhoTheta_comp, 0, 1, ng);
+
1038  MultiFab::Divide(*Theta_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);
+
1039 
+
1040  if (solverChoice.moisture_type != MoistureType::None) {
+
1041  ng = Qv_prim[lev]->nGrowVect();
+
1042 
+
1043  MultiFab::Copy( *Qv_prim[lev], vars_new[lev][Vars::cons], RhoQ1_comp, 0, 1, ng);
+
1044  MultiFab::Divide(*Qv_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);
+
1045 
+
1046  int rhoqr_comp = solverChoice.RhoQr_comp;
+
1047  if (rhoqr_comp > -1) {
+
1048  MultiFab::Copy( *Qr_prim[lev], vars_new[lev][Vars::cons], rhoqr_comp, 0, 1, ng);
+
1049  MultiFab::Divide(*Qr_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);
+
1050  } else {
+
1051  Qr_prim[lev]->setVal(0.0);
+
1052  }
+
1053  }
+
1054  m_most->update_mac_ptrs(lev, vars_new, Theta_prim, Qv_prim, Qr_prim);
+
1055 
+
1056  if (restart_chkfile == "") {
+
1057  // Only do this if starting from scratch; if restarting, then
+
1058  // we don't want to call update_fluxes multiple times because
+
1059  // it will change u* and theta* from their previous values
+
1060  m_most->update_pblh(lev, vars_new, z_phys_cc[lev].get(),
+ + + +
1064  m_most->update_fluxes(lev, time);
+
1065  }
+
1066  }
+
1067  }
+
1068 
+
1069  // Update micro vars before first plot file
+
1070  if (solverChoice.moisture_type != MoistureType::None) {
+
1071  for (int lev = 0; lev <= finest_level; ++lev) micro->Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
+
1072  }
+
1073 
+
1074  // Fill time averaged velocities before first plot file
+
1075  if (solverChoice.time_avg_vel) {
+
1076  for (int lev = 0; lev <= finest_level; ++lev) {
+
1077  Time_Avg_Vel_atCC(dt[lev], t_avg_cnt[lev], vel_t_avg[lev].get(),
+
1078  vars_new[lev][Vars::xvel],
+
1079  vars_new[lev][Vars::yvel],
+
1080  vars_new[lev][Vars::zvel]);
1081  }
- -
1083  }
-
1084 
-
1085  if ( (restart_chkfile.empty()) ||
-
1086  (!restart_chkfile.empty() && plot_file_on_restart) )
-
1087  {
-
1088  if (m_plot_int_1 > 0 || m_plot_per_1 > 0.)
-
1089  {
- - -
1092  }
-
1093  if (m_plot_int_2 > 0 || m_plot_per_2 > 0.)
-
1094  {
- - -
1097  }
-
1098  }
-
1099 
-
1100  // Set these up here because we need to know which MPI rank "cell" is on...
-
1101  if (pp.contains("data_log"))
-
1102  {
-
1103  int num_datalogs = pp.countval("data_log");
-
1104  datalog.resize(num_datalogs);
-
1105  datalogname.resize(num_datalogs);
-
1106  pp.queryarr("data_log",datalogname,0,num_datalogs);
-
1107  for (int i = 0; i < num_datalogs; i++)
- -
1109  }
-
1110 
-
1111  if (restart_chkfile.empty() && profile_int > 0) {
-
1112  if (destag_profiles) {
-
1113  // all variables cell-centered
- -
1115  } else {
-
1116  // some variables staggered
- -
1118  }
-
1119  }
-
1120 
-
1121  if (pp.contains("sample_point_log") && pp.contains("sample_point"))
-
1122  {
-
1123  int lev = 0;
-
1124 
-
1125  int num_samplepts = pp.countval("sample_point") / AMREX_SPACEDIM;
-
1126  if (num_samplepts > 0) {
-
1127  Vector<int> index; index.resize(num_samplepts*AMREX_SPACEDIM);
-
1128  samplepoint.resize(num_samplepts);
-
1129 
-
1130  pp.queryarr("sample_point",index,0,num_samplepts*AMREX_SPACEDIM);
-
1131  for (int i = 0; i < num_samplepts; i++) {
-
1132  IntVect iv(index[AMREX_SPACEDIM*i+0],index[AMREX_SPACEDIM*i+1],index[AMREX_SPACEDIM*i+2]);
-
1133  samplepoint[i] = iv;
-
1134  }
+
1082  }
+
1083 
+
1084  // check for additional plotting variables that are available after particle containers
+
1085  // are setup.
+
1086  const std::string& pv1 = "plot_vars_1"; appendPlotVariables(pv1,plot_var_names_1);
+
1087  const std::string& pv2 = "plot_vars_2"; appendPlotVariables(pv2,plot_var_names_2);
+
1088 
+
1089  if ( restart_chkfile.empty() && (m_check_int > 0 || m_check_per > 0.) )
+
1090  {
+
1091 #ifdef ERF_USE_NETCDF
+
1092  if (check_type == "netcdf") {
+
1093  WriteNCCheckpointFile();
+
1094  }
+
1095 #endif
+
1096  if (check_type == "native") {
+ +
1098  }
+ +
1100  }
+
1101 
+
1102  if ( (restart_chkfile.empty()) ||
+
1103  (!restart_chkfile.empty() && plot_file_on_restart) )
+
1104  {
+
1105  if (m_plot_int_1 > 0 || m_plot_per_1 > 0.)
+
1106  {
+ + +
1109  }
+
1110  if (m_plot_int_2 > 0 || m_plot_per_2 > 0.)
+
1111  {
+ + +
1114  }
+
1115  }
+
1116 
+
1117  // Set these up here because we need to know which MPI rank "cell" is on...
+
1118  if (pp.contains("data_log"))
+
1119  {
+
1120  int num_datalogs = pp.countval("data_log");
+
1121  datalog.resize(num_datalogs);
+
1122  datalogname.resize(num_datalogs);
+
1123  pp.queryarr("data_log",datalogname,0,num_datalogs);
+
1124  for (int i = 0; i < num_datalogs; i++)
+ +
1126  }
+
1127 
+
1128  if (restart_chkfile.empty() && profile_int > 0) {
+
1129  if (destag_profiles) {
+
1130  // all variables cell-centered
+ +
1132  } else {
+
1133  // some variables staggered
+
1135  }
-
1136 
-
1137  int num_sampleptlogs = pp.countval("sample_point_log");
-
1138  AMREX_ALWAYS_ASSERT(num_sampleptlogs == num_samplepts);
-
1139  if (num_sampleptlogs > 0) {
-
1140  sampleptlog.resize(num_sampleptlogs);
-
1141  sampleptlogname.resize(num_sampleptlogs);
-
1142  pp.queryarr("sample_point_log",sampleptlogname,0,num_sampleptlogs);
-
1143 
-
1144  for (int i = 0; i < num_sampleptlogs; i++) {
- -
1146  }
-
1147  }
-
1148 
-
1149  }
-
1150 
-
1151  if (pp.contains("sample_line_log") && pp.contains("sample_line"))
-
1152  {
-
1153  int lev = 0;
-
1154 
-
1155  int num_samplelines = pp.countval("sample_line") / AMREX_SPACEDIM;
-
1156  if (num_samplelines > 0) {
-
1157  Vector<int> index; index.resize(num_samplelines*AMREX_SPACEDIM);
-
1158  sampleline.resize(num_samplelines);
-
1159 
-
1160  pp.queryarr("sample_line",index,0,num_samplelines*AMREX_SPACEDIM);
-
1161  for (int i = 0; i < num_samplelines; i++) {
-
1162  IntVect iv(index[AMREX_SPACEDIM*i+0],index[AMREX_SPACEDIM*i+1],index[AMREX_SPACEDIM*i+2]);
-
1163  sampleline[i] = iv;
-
1164  }
-
1165  }
-
1166 
-
1167  int num_samplelinelogs = pp.countval("sample_line_log");
-
1168  AMREX_ALWAYS_ASSERT(num_samplelinelogs == num_samplelines);
-
1169  if (num_samplelinelogs > 0) {
-
1170  samplelinelog.resize(num_samplelinelogs);
-
1171  samplelinelogname.resize(num_samplelinelogs);
-
1172  pp.queryarr("sample_line_log",samplelinelogname,0,num_samplelinelogs);
-
1173 
-
1174  for (int i = 0; i < num_samplelinelogs; i++) {
- -
1176  }
-
1177  }
-
1178 
-
1179  }
-
1180 
-
1181  // Create object to do line and plane sampling if needed
-
1182  bool do_line = false; bool do_plane = false;
-
1183  pp.query("do_line_sampling",do_line); pp.query("do_plane_sampling",do_plane);
-
1184  if (do_line || do_plane) { data_sampler = std::make_unique<SampleData>(do_line, do_plane); }
-
1185 
-
1186  if ( solverChoice.terrain_type == TerrainType::EB ||
-
1187  solverChoice.terrain_type == TerrainType::ImmersedForcing)
-
1188  {
-
1189  bool write_eb_surface = false;
-
1190  pp.query("write_eb_surface", write_eb_surface);
-
1191  if (write_eb_surface) WriteMyEBSurface();
-
1192  }
-
1193 
-
1194 }
+
1136  }
+
1137 
+
1138  if (pp.contains("sample_point_log") && pp.contains("sample_point"))
+
1139  {
+
1140  int lev = 0;
+
1141 
+
1142  int num_samplepts = pp.countval("sample_point") / AMREX_SPACEDIM;
+
1143  if (num_samplepts > 0) {
+
1144  Vector<int> index; index.resize(num_samplepts*AMREX_SPACEDIM);
+
1145  samplepoint.resize(num_samplepts);
+
1146 
+
1147  pp.queryarr("sample_point",index,0,num_samplepts*AMREX_SPACEDIM);
+
1148  for (int i = 0; i < num_samplepts; i++) {
+
1149  IntVect iv(index[AMREX_SPACEDIM*i+0],index[AMREX_SPACEDIM*i+1],index[AMREX_SPACEDIM*i+2]);
+
1150  samplepoint[i] = iv;
+
1151  }
+
1152  }
+
1153 
+
1154  int num_sampleptlogs = pp.countval("sample_point_log");
+
1155  AMREX_ALWAYS_ASSERT(num_sampleptlogs == num_samplepts);
+
1156  if (num_sampleptlogs > 0) {
+
1157  sampleptlog.resize(num_sampleptlogs);
+
1158  sampleptlogname.resize(num_sampleptlogs);
+
1159  pp.queryarr("sample_point_log",sampleptlogname,0,num_sampleptlogs);
+
1160 
+
1161  for (int i = 0; i < num_sampleptlogs; i++) {
+ +
1163  }
+
1164  }
+
1165 
+
1166  }
+
1167 
+
1168  if (pp.contains("sample_line_log") && pp.contains("sample_line"))
+
1169  {
+
1170  int lev = 0;
+
1171 
+
1172  int num_samplelines = pp.countval("sample_line") / AMREX_SPACEDIM;
+
1173  if (num_samplelines > 0) {
+
1174  Vector<int> index; index.resize(num_samplelines*AMREX_SPACEDIM);
+
1175  sampleline.resize(num_samplelines);
+
1176 
+
1177  pp.queryarr("sample_line",index,0,num_samplelines*AMREX_SPACEDIM);
+
1178  for (int i = 0; i < num_samplelines; i++) {
+
1179  IntVect iv(index[AMREX_SPACEDIM*i+0],index[AMREX_SPACEDIM*i+1],index[AMREX_SPACEDIM*i+2]);
+
1180  sampleline[i] = iv;
+
1181  }
+
1182  }
+
1183 
+
1184  int num_samplelinelogs = pp.countval("sample_line_log");
+
1185  AMREX_ALWAYS_ASSERT(num_samplelinelogs == num_samplelines);
+
1186  if (num_samplelinelogs > 0) {
+
1187  samplelinelog.resize(num_samplelinelogs);
+
1188  samplelinelogname.resize(num_samplelinelogs);
+
1189  pp.queryarr("sample_line_log",samplelinelogname,0,num_samplelinelogs);
+
1190 
+
1191  for (int i = 0; i < num_samplelinelogs; i++) {
+ +
1193  }
+
1194  }
+
1195 
+
1196  }
+
1197 
+
1198  // Create object to do line and plane sampling if needed
+
1199  bool do_line = false; bool do_plane = false;
+
1200  pp.query("do_line_sampling",do_line); pp.query("do_plane_sampling",do_plane);
+
1201  if (do_line || do_plane) { data_sampler = std::make_unique<SampleData>(do_line, do_plane); }
+
1202 
+
1203  if ( solverChoice.terrain_type == TerrainType::EB ||
+
1204  solverChoice.terrain_type == TerrainType::ImmersedForcing)
+
1205  {
+
1206  bool write_eb_surface = false;
+
1207  pp.query("write_eb_surface", write_eb_surface);
+
1208  if (write_eb_surface) WriteMyEBSurface();
+
1209  }
+
1210 
+
1211 }
void initRayleigh()
Initialize Rayleigh damping profiles.
Definition: ERF_InitRayleigh.cpp:14
amrex::Vector< std::string > samplelinelogname
Definition: ERF.H:1382
void setRayleighRefFromSounding(bool restarting)
Set Rayleigh mean profiles from input sounding.
Definition: ERF_InitRayleigh.cpp:55
@@ -9587,7 +9602,7 @@

amrex::Vector< std::unique_ptr< std::fstream > > samplelinelog
Definition: ERF.H:1381
static int sum_interval
Definition: ERF.H:1047
static int pert_interval
Definition: ERF.H:1048
-
void restart()
Definition: ERF.cpp:1234
+
void restart()
Definition: ERF.cpp:1251
void ReadCheckpointFileMOST()
Definition: ERF_Checkpoint.cpp:627
void write_1D_profiles(amrex::Real time)
Definition: ERF_Write1DProfiles.cpp:15
int profile_int
Definition: ERF.H:952
@@ -9607,7 +9622,7 @@

static bool is_it_time_for_action(int nstep, amrex::Real time, amrex::Real dt, int action_interval, amrex::Real action_per)
Definition: ERF_WriteScalarProfiles.cpp:466
void project_velocities(int lev, amrex::Real dt, amrex::Vector< amrex::MultiFab > &vars, amrex::MultiFab &p)
Definition: ERF_PoissonSolve.cpp:10
int plot_file_on_restart
Definition: ERF.H:911
-
void Construct_ERFFillPatchers(int lev)
Definition: ERF.cpp:1876
+
void Construct_ERFFillPatchers(int lev)
Definition: ERF.cpp:1893
void setRecordSampleLineInfo(int i, int lev, amrex::IntVect &cell, const std::string &filename)
Definition: ERF.H:1352
void setSpongeRefFromSounding(bool restarting)
Set sponge mean profiles from input sounding.
Definition: ERF_InitSponge.cpp:65
amrex::Vector< amrex::IntVect > samplepoint
Definition: ERF.H:1379
@@ -9913,25 +9928,25 @@

-
591 {
-
592  const BoxArray& ba(cons_mf.boxArray());
-
593  const DistributionMapping& dm(cons_mf.DistributionMap());
+
589 {
+
590  const BoxArray& ba(cons_mf.boxArray());
+
591  const DistributionMapping& dm(cons_mf.DistributionMap());
+
592 
+
593  int ncomp_cons = cons_mf.nComp();
594 
-
595  int ncomp_cons = cons_mf.nComp();
-
596 
-
597  // Initialize the integrator memory
-
598  Vector<MultiFab> int_state; // integration state data structure example
-
599  int_state.push_back(MultiFab(cons_mf, make_alias, 0, ncomp_cons)); // cons
-
600  int_state.push_back(MultiFab(convert(ba,IntVect(1,0,0)), dm, 1, vel_mf.nGrow())); // xmom
-
601  int_state.push_back(MultiFab(convert(ba,IntVect(0,1,0)), dm, 1, vel_mf.nGrow())); // ymom
-
602  int_state.push_back(MultiFab(convert(ba,IntVect(0,0,1)), dm, 1, vel_mf.nGrow())); // zmom
-
603 
-
604  mri_integrator_mem[lev] = std::make_unique<MRISplitIntegrator<Vector<MultiFab> > >(int_state);
-
605  mri_integrator_mem[lev]->setNoSubstepping((solverChoice.substepping_type[lev] == SubsteppingType::None));
-
606  mri_integrator_mem[lev]->setAnelastic(solverChoice.anelastic[lev]);
-
607  mri_integrator_mem[lev]->setNcompCons(ncomp_cons);
-
608  mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(solverChoice.force_stage1_single_substep);
-
609 }
+
595  // Initialize the integrator memory
+
596  Vector<MultiFab> int_state; // integration state data structure example
+
597  int_state.push_back(MultiFab(cons_mf, make_alias, 0, ncomp_cons)); // cons
+
598  int_state.push_back(MultiFab(convert(ba,IntVect(1,0,0)), dm, 1, vel_mf.nGrow())); // xmom
+
599  int_state.push_back(MultiFab(convert(ba,IntVect(0,1,0)), dm, 1, vel_mf.nGrow())); // ymom
+
600  int_state.push_back(MultiFab(convert(ba,IntVect(0,0,1)), dm, 1, vel_mf.nGrow())); // zmom
+
601 
+
602  mri_integrator_mem[lev] = std::make_unique<MRISplitIntegrator<Vector<MultiFab> > >(int_state);
+
603  mri_integrator_mem[lev]->setNoSubstepping((solverChoice.substepping_type[lev] == SubsteppingType::None));
+
604  mri_integrator_mem[lev]->setAnelastic(solverChoice.anelastic[lev]);
+
605  mri_integrator_mem[lev]->setNcompCons(ncomp_cons);
+
606  mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(solverChoice.force_stage1_single_substep);
+
607 }
@@ -10025,30 +10040,30 @@

1199 {
-
1200  if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Eulerian) {
-
1201 
-
1202  micro = std::make_unique<EulerianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
-
1203 
-
1204  } else if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
-
1205 #ifdef ERF_USE_PARTICLES
-
1206 
-
1207  micro = std::make_unique<LagrangianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
-
1208  /* Lagrangian microphysics models will have a particle container; it needs to be added
-
1209  to ERF::particleData */
-
1210  const auto& pc_name( dynamic_cast<LagrangianMicrophysics&>(*micro).getName() );
-
1211  /* The particle container has not yet been constructed and initialized, so just add
-
1212  its name here for now (so that functions to set plotting variables can see it). */
-
1213  particleData.addName( pc_name );
-
1214 
-
1215 #else
-
1216  Abort("Lagrangian microphysics can be used when compiled with ERF_USE_PARTICLES");
-
1217 #endif
-
1218  }
-
1219 
-
1220  qmoist.resize(a_nlevsmax);
-
1221  return;
-
1222 }
+
1216 {
+
1217  if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Eulerian) {
+
1218 
+
1219  micro = std::make_unique<EulerianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
+
1220 
+
1221  } else if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
+
1222 #ifdef ERF_USE_PARTICLES
+
1223 
+
1224  micro = std::make_unique<LagrangianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
+
1225  /* Lagrangian microphysics models will have a particle container; it needs to be added
+
1226  to ERF::particleData */
+
1227  const auto& pc_name( dynamic_cast<LagrangianMicrophysics&>(*micro).getName() );
+
1228  /* The particle container has not yet been constructed and initialized, so just add
+
1229  its name here for now (so that functions to set plotting variables can see it). */
+
1230  particleData.addName( pc_name );
+
1231 
+
1232 #else
+
1233  Abort("Lagrangian microphysics can be used when compiled with ERF_USE_PARTICLES");
+
1234 #endif
+
1235  }
+
1236 
+
1237  qmoist.resize(a_nlevsmax);
+
1238  return;
+
1239 }
amrex::Vector< amrex::Vector< amrex::MultiFab * > > qmoist
Definition: ERF.H:769
Here is the call graph for this function:
@@ -10489,27 +10504,27 @@

-
613 {
-
614  if (SolverChoice::mesh_type == MeshType::VariableDz) {
-
615  AMREX_ALWAYS_ASSERT(z_phys_nd[lev] != nullptr);
-
616  }
-
617 
-
618  physbcs_cons[lev] = std::make_unique<ERFPhysBCFunct_cons> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
- -
620  z_phys_nd[lev], use_real_bcs);
-
621  physbcs_u[lev] = std::make_unique<ERFPhysBCFunct_u> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
- -
623  z_phys_nd[lev], use_real_bcs, xvel_bc_data[lev].data());
-
624  physbcs_v[lev] = std::make_unique<ERFPhysBCFunct_v> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
- -
626  z_phys_nd[lev], use_real_bcs, yvel_bc_data[lev].data());
-
627  physbcs_w[lev] = std::make_unique<ERFPhysBCFunct_w> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
- - -
630  use_real_bcs, zvel_bc_data[lev].data());
-
631  physbcs_base[lev] = std::make_unique<ERFPhysBCFunct_base> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
-
632  (solverChoice.terrain_type == TerrainType::MovingFittedMesh));
-
633 }
+
611 {
+
612  if (SolverChoice::mesh_type == MeshType::VariableDz) {
+
613  AMREX_ALWAYS_ASSERT(z_phys_nd[lev] != nullptr);
+
614  }
+
615 
+
616  physbcs_cons[lev] = std::make_unique<ERFPhysBCFunct_cons> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
+ +
618  z_phys_nd[lev], use_real_bcs);
+
619  physbcs_u[lev] = std::make_unique<ERFPhysBCFunct_u> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
+ +
621  z_phys_nd[lev], use_real_bcs, xvel_bc_data[lev].data());
+
622  physbcs_v[lev] = std::make_unique<ERFPhysBCFunct_v> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
+ +
624  z_phys_nd[lev], use_real_bcs, yvel_bc_data[lev].data());
+
625  physbcs_w[lev] = std::make_unique<ERFPhysBCFunct_w> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
+ + +
628  use_real_bcs, zvel_bc_data[lev].data());
+
629  physbcs_base[lev] = std::make_unique<ERFPhysBCFunct_base> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
+
630  (solverChoice.terrain_type == TerrainType::MovingFittedMesh));
+
631 }
@@ -10544,49 +10559,49 @@

-
1831 {
-
1832  // Get the number of cells in z at level 0
-
1833  int dir_z = AMREX_SPACEDIM-1;
-
1834  auto domain = geom[0].Domain();
-
1835  int size_z = domain.length(dir_z);
-
1836  int start_z = domain.smallEnd()[dir_z];
-
1837  Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
-
1838 
-
1839  // resize the level 0 horizontal average vectors
-
1840  h_havg.resize(size_z, 0.0_rt);
-
1841 
-
1842  // Get the cell centered data and construct sums
-
1843 #ifdef _OPENMP
-
1844 #pragma omp parallel if (Gpu::notInLaunchRegion())
-
1845 #endif
-
1846  for (MFIter mfi(S); mfi.isValid(); ++mfi) {
-
1847  const Box& box = mfi.validbox();
-
1848  const IntVect& se = box.smallEnd();
-
1849  const IntVect& be = box.bigEnd();
-
1850 
-
1851  auto fab_arr = S[mfi].array();
-
1852 
-
1853  FArrayBox fab_reduce(box, 1, The_Async_Arena());
-
1854  auto arr_reduce = fab_reduce.array();
+
1848 {
+
1849  // Get the number of cells in z at level 0
+
1850  int dir_z = AMREX_SPACEDIM-1;
+
1851  auto domain = geom[0].Domain();
+
1852  int size_z = domain.length(dir_z);
+
1853  int start_z = domain.smallEnd()[dir_z];
+
1854  Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
1855 
-
1856  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-
1857  arr_reduce(i, j, k, 0) = fab_arr(i,j,k,n);
-
1858  });
-
1859 
-
1860  for (int k=se[dir_z]; k <= be[dir_z]; ++k) {
-
1861  Box kbox(box); kbox.setSmall(dir_z,k); kbox.setBig(dir_z,k);
-
1862  h_havg[k-start_z] += fab_reduce.sum<RunOn::Device>(kbox,0);
-
1863  }
-
1864  }
-
1865 
-
1866  // combine sums from different MPI ranks
-
1867  ParallelDescriptor::ReduceRealSum(h_havg.dataPtr(), h_havg.size());
-
1868 
-
1869  // divide by the total number of cells we are averaging over
-
1870  for (int k = 0; k < size_z; ++k) {
-
1871  h_havg[k] /= area_z;
-
1872  }
-
1873 }
+
1856  // resize the level 0 horizontal average vectors
+
1857  h_havg.resize(size_z, 0.0_rt);
+
1858 
+
1859  // Get the cell centered data and construct sums
+
1860 #ifdef _OPENMP
+
1861 #pragma omp parallel if (Gpu::notInLaunchRegion())
+
1862 #endif
+
1863  for (MFIter mfi(S); mfi.isValid(); ++mfi) {
+
1864  const Box& box = mfi.validbox();
+
1865  const IntVect& se = box.smallEnd();
+
1866  const IntVect& be = box.bigEnd();
+
1867 
+
1868  auto fab_arr = S[mfi].array();
+
1869 
+
1870  FArrayBox fab_reduce(box, 1, The_Async_Arena());
+
1871  auto arr_reduce = fab_reduce.array();
+
1872 
+
1873  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+
1874  arr_reduce(i, j, k, 0) = fab_arr(i,j,k,n);
+
1875  });
+
1876 
+
1877  for (int k=se[dir_z]; k <= be[dir_z]; ++k) {
+
1878  Box kbox(box); kbox.setSmall(dir_z,k); kbox.setBig(dir_z,k);
+
1879  h_havg[k-start_z] += fab_reduce.sum<RunOn::Device>(kbox,0);
+
1880  }
+
1881  }
+
1882 
+
1883  // combine sums from different MPI ranks
+
1884  ParallelDescriptor::ReduceRealSum(h_havg.dataPtr(), h_havg.size());
+
1885 
+
1886  // divide by the total number of cells we are averaging over
+
1887  for (int k = 0; k < size_z; ++k) {
+
1888  h_havg[k] /= area_z;
+
1889  }
+
1890 }
@@ -10684,106 +10699,106 @@

-
1725 {
-
1726  int lev = 0;
-
1727 
-
1728  // First, average down all levels (if doing two-way coupling)
-
1729  if (solverChoice.coupling_type == CouplingType::TwoWay) {
-
1730  AverageDown();
-
1731  }
-
1732 
-
1733  MultiFab mf(grids[lev], dmap[lev], 5, 0);
-
1734 
-
1735  int zdir = 2;
-
1736  auto domain = geom[0].Domain();
-
1737 
-
1738  bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
-
1739  bool is_anelastic = (solverChoice.anelastic[lev] == 1);
-
1740 
-
1741  for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
-
1742  const Box& bx = mfi.validbox();
-
1743  auto fab_arr = mf.array(mfi);
-
1744  auto const hse_arr = base_state[lev].const_array(mfi);
-
1745  auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
-
1746  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-
1747  Real dens = cons_arr(i, j, k, Rho_comp);
-
1748  fab_arr(i, j, k, 0) = dens;
-
1749  fab_arr(i, j, k, 1) = cons_arr(i, j, k, RhoTheta_comp) / dens;
-
1750  if (!use_moisture) {
-
1751  if (is_anelastic) {
-
1752  fab_arr(i,j,k,2) = hse_arr(i,j,k,BaseState::p0_comp);
-
1753  } else {
-
1754  fab_arr(i,j,k,2) = getPgivenRTh(cons_arr(i,j,k,RhoTheta_comp));
-
1755  }
-
1756  }
-
1757  });
-
1758  }
-
1759 
-
1760  if (use_moisture)
-
1761  {
-
1762  for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
-
1763  const Box& bx = mfi.validbox();
-
1764  auto fab_arr = mf.array(mfi);
-
1765  auto const hse_arr = base_state[lev].const_array(mfi);
-
1766  auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
-
1767  int ncomp = vars_new[lev][Vars::cons].nComp();
-
1768 
-
1769  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-
1770  Real dens = cons_arr(i, j, k, Rho_comp);
-
1771  if (is_anelastic) {
-
1772  fab_arr(i,j,k,2) = hse_arr(i,j,k,BaseState::p0_comp);
-
1773  } else {
-
1774  Real qv = cons_arr(i, j, k, RhoQ1_comp) / dens;
-
1775  fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv);
-
1776  }
-
1777  fab_arr(i, j, k, 3) = (ncomp > RhoQ1_comp ? cons_arr(i, j, k, RhoQ1_comp) / dens : 0.0);
-
1778  fab_arr(i, j, k, 4) = (ncomp > RhoQ2_comp ? cons_arr(i, j, k, RhoQ2_comp) / dens : 0.0);
-
1779  });
-
1780  }
-
1781 
-
1782  Gpu::HostVector<Real> h_avg_qv = sumToLine(mf,3,1,domain,zdir);
-
1783  Gpu::HostVector<Real> h_avg_qc = sumToLine(mf,4,1,domain,zdir);
-
1784  }
+
1742 {
+
1743  int lev = 0;
+
1744 
+
1745  // First, average down all levels (if doing two-way coupling)
+
1746  if (solverChoice.coupling_type == CouplingType::TwoWay) {
+
1747  AverageDown();
+
1748  }
+
1749 
+
1750  MultiFab mf(grids[lev], dmap[lev], 5, 0);
+
1751 
+
1752  int zdir = 2;
+
1753  auto domain = geom[0].Domain();
+
1754 
+
1755  bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
+
1756  bool is_anelastic = (solverChoice.anelastic[lev] == 1);
+
1757 
+
1758  for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
+
1759  const Box& bx = mfi.validbox();
+
1760  auto fab_arr = mf.array(mfi);
+
1761  auto const hse_arr = base_state[lev].const_array(mfi);
+
1762  auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
+
1763  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+
1764  Real dens = cons_arr(i, j, k, Rho_comp);
+
1765  fab_arr(i, j, k, 0) = dens;
+
1766  fab_arr(i, j, k, 1) = cons_arr(i, j, k, RhoTheta_comp) / dens;
+
1767  if (!use_moisture) {
+
1768  if (is_anelastic) {
+
1769  fab_arr(i,j,k,2) = hse_arr(i,j,k,BaseState::p0_comp);
+
1770  } else {
+
1771  fab_arr(i,j,k,2) = getPgivenRTh(cons_arr(i,j,k,RhoTheta_comp));
+
1772  }
+
1773  }
+
1774  });
+
1775  }
+
1776 
+
1777  if (use_moisture)
+
1778  {
+
1779  for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
+
1780  const Box& bx = mfi.validbox();
+
1781  auto fab_arr = mf.array(mfi);
+
1782  auto const hse_arr = base_state[lev].const_array(mfi);
+
1783  auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
+
1784  int ncomp = vars_new[lev][Vars::cons].nComp();
1785 
-
1786  // Sum in the horizontal plane
-
1787  Gpu::HostVector<Real> h_avg_density = sumToLine(mf,0,1,domain,zdir);
-
1788  Gpu::HostVector<Real> h_avg_temperature = sumToLine(mf,1,1,domain,zdir);
-
1789  Gpu::HostVector<Real> h_avg_pressure = sumToLine(mf,2,1,domain,zdir);
-
1790 
-
1791  // Divide by the total number of cells we are averaging over
-
1792  int size_z = domain.length(zdir);
-
1793  Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
-
1794  int klen = static_cast<int>(h_avg_density.size());
-
1795 
-
1796  for (int k = 0; k < klen; ++k) {
-
1797  h_havg_density[k] /= area_z;
-
1798  h_havg_temperature[k] /= area_z;
-
1799  h_havg_pressure[k] /= area_z;
-
1800  if (solverChoice.moisture_type != MoistureType::None)
-
1801  {
-
1802  h_havg_qc[k] /= area_z;
-
1803  h_havg_qv[k] /= area_z;
-
1804  }
-
1805  } // k
-
1806 
-
1807  // resize device vectors
-
1808  d_havg_density.resize(size_z, 0.0_rt);
-
1809  d_havg_temperature.resize(size_z, 0.0_rt);
-
1810  d_havg_pressure.resize(size_z, 0.0_rt);
-
1811 
-
1812  // copy host vectors to device vectors
-
1813  Gpu::copy(Gpu::hostToDevice, h_havg_density.begin(), h_havg_density.end(), d_havg_density.begin());
-
1814  Gpu::copy(Gpu::hostToDevice, h_havg_temperature.begin(), h_havg_temperature.end(), d_havg_temperature.begin());
-
1815  Gpu::copy(Gpu::hostToDevice, h_havg_pressure.begin(), h_havg_pressure.end(), d_havg_pressure.begin());
-
1816 
-
1817  if (solverChoice.moisture_type != MoistureType::None)
-
1818  {
-
1819  d_havg_qv.resize(size_z, 0.0_rt);
-
1820  d_havg_qc.resize(size_z, 0.0_rt);
-
1821  Gpu::copy(Gpu::hostToDevice, h_havg_qv.begin(), h_havg_qv.end(), d_havg_qv.begin());
-
1822  Gpu::copy(Gpu::hostToDevice, h_havg_qc.begin(), h_havg_qc.end(), d_havg_qc.begin());
-
1823  }
-
1824 }
+
1786  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+
1787  Real dens = cons_arr(i, j, k, Rho_comp);
+
1788  if (is_anelastic) {
+
1789  fab_arr(i,j,k,2) = hse_arr(i,j,k,BaseState::p0_comp);
+
1790  } else {
+
1791  Real qv = cons_arr(i, j, k, RhoQ1_comp) / dens;
+
1792  fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv);
+
1793  }
+
1794  fab_arr(i, j, k, 3) = (ncomp > RhoQ1_comp ? cons_arr(i, j, k, RhoQ1_comp) / dens : 0.0);
+
1795  fab_arr(i, j, k, 4) = (ncomp > RhoQ2_comp ? cons_arr(i, j, k, RhoQ2_comp) / dens : 0.0);
+
1796  });
+
1797  }
+
1798 
+
1799  Gpu::HostVector<Real> h_avg_qv = sumToLine(mf,3,1,domain,zdir);
+
1800  Gpu::HostVector<Real> h_avg_qc = sumToLine(mf,4,1,domain,zdir);
+
1801  }
+
1802 
+
1803  // Sum in the horizontal plane
+
1804  Gpu::HostVector<Real> h_avg_density = sumToLine(mf,0,1,domain,zdir);
+
1805  Gpu::HostVector<Real> h_avg_temperature = sumToLine(mf,1,1,domain,zdir);
+
1806  Gpu::HostVector<Real> h_avg_pressure = sumToLine(mf,2,1,domain,zdir);
+
1807 
+
1808  // Divide by the total number of cells we are averaging over
+
1809  int size_z = domain.length(zdir);
+
1810  Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
+
1811  int klen = static_cast<int>(h_avg_density.size());
+
1812 
+
1813  for (int k = 0; k < klen; ++k) {
+
1814  h_havg_density[k] /= area_z;
+
1815  h_havg_temperature[k] /= area_z;
+
1816  h_havg_pressure[k] /= area_z;
+
1817  if (solverChoice.moisture_type != MoistureType::None)
+
1818  {
+
1819  h_havg_qc[k] /= area_z;
+
1820  h_havg_qv[k] /= area_z;
+
1821  }
+
1822  } // k
+
1823 
+
1824  // resize device vectors
+
1825  d_havg_density.resize(size_z, 0.0_rt);
+
1826  d_havg_temperature.resize(size_z, 0.0_rt);
+
1827  d_havg_pressure.resize(size_z, 0.0_rt);
+
1828 
+
1829  // copy host vectors to device vectors
+
1830  Gpu::copy(Gpu::hostToDevice, h_havg_density.begin(), h_havg_density.end(), d_havg_density.begin());
+
1831  Gpu::copy(Gpu::hostToDevice, h_havg_temperature.begin(), h_havg_temperature.end(), d_havg_temperature.begin());
+
1832  Gpu::copy(Gpu::hostToDevice, h_havg_pressure.begin(), h_havg_pressure.end(), d_havg_pressure.begin());
+
1833 
+
1834  if (solverChoice.moisture_type != MoistureType::None)
+
1835  {
+
1836  d_havg_qv.resize(size_z, 0.0_rt);
+
1837  d_havg_qc.resize(size_z, 0.0_rt);
+
1838  Gpu::copy(Gpu::hostToDevice, h_havg_qv.begin(), h_havg_qv.end(), d_havg_qv.begin());
+
1839  Gpu::copy(Gpu::hostToDevice, h_havg_qc.begin(), h_havg_qc.end(), d_havg_qc.begin());
+
1840  }
+
1841 }
amrex::Gpu::DeviceVector< amrex::Real > d_havg_temperature
Definition: ERF.H:1156
amrex::Gpu::DeviceVector< amrex::Real > d_havg_qv
Definition: ERF.H:1158
amrex::Vector< amrex::Real > h_havg_pressure
Definition: ERF.H:1151
@@ -10959,12 +10974,12 @@

316  // particleData.Redistribute();

317 #endif
318 }
-
void update_diffusive_arrays(int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
Definition: ERF_MakeNewArrays.cpp:382
-
void initialize_integrator(int lev, amrex::MultiFab &cons_mf, amrex::MultiFab &vel_mf)
Definition: ERF_MakeNewArrays.cpp:590
-
void update_terrain_arrays(int lev)
Definition: ERF_MakeNewArrays.cpp:579
-
void init_zphys(int lev, amrex::Real time)
Definition: ERF_MakeNewArrays.cpp:473
+
void update_diffusive_arrays(int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
Definition: ERF_MakeNewArrays.cpp:380
+
void initialize_integrator(int lev, amrex::MultiFab &cons_mf, amrex::MultiFab &vel_mf)
Definition: ERF_MakeNewArrays.cpp:588
+
void update_terrain_arrays(int lev)
Definition: ERF_MakeNewArrays.cpp:577
+
void init_zphys(int lev, amrex::Real time)
Definition: ERF_MakeNewArrays.cpp:471
void init_stuff(int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm, amrex::Vector< amrex::MultiFab > &lev_new, amrex::Vector< amrex::MultiFab > &lev_old, amrex::MultiFab &tmp_base_state, std::unique_ptr< amrex::MultiFab > &tmp_zphys_nd)
Definition: ERF_MakeNewArrays.cpp:23
-
void Define_ERFFillPatchers(int lev)
Definition: ERF.cpp:1902
+
void Define_ERFFillPatchers(int lev)
Definition: ERF.cpp:1919
bool do_forest_drag
Definition: ERF_DataStruct.H:710
@@ -11195,7 +11210,7 @@

BoxArray ERFPostProcessBaseGrids(const Box &domain, bool decompose_in_z)
Definition: ERF_ChopGrids.cpp:6
void thinbody_wall_dist(std::unique_ptr< MultiFab > &wdist, Vector< IntVect > &xfaces, Vector< IntVect > &yfaces, Vector< IntVect > &zfaces, const Geometry &geomdata, std::unique_ptr< MultiFab > &z_phys_cc)
Definition: ERF_ThinBodyWallDist.cpp:12
static int nghost_eb_basic()
Definition: ERF.H:1404
-
void init_only(int lev, amrex::Real time)
Definition: ERF.cpp:1265
+
void init_only(int lev, amrex::Real time)
Definition: ERF.cpp:1282
void poisson_wall_dist(int lev)
Definition: ERF_PoissonWallDist.cpp:20
static int nghost_eb_volume()
Definition: ERF.H:1408
static int nghost_eb_full()
Definition: ERF.H:1411
@@ -11510,67 +11525,67 @@

-
1659 {
-
1660  AMREX_ALWAYS_ASSERT(cfl > 0. || fixed_dt[0] > 0.);
-
1661 
-
1662  // We don't allow use_real_bcs to be true if init_type is not either InitType::Real or InitType::Metgrid
-
1663  AMREX_ALWAYS_ASSERT(!use_real_bcs || ((init_type == InitType::Real) || (init_type == InitType::Metgrid)) );
-
1664 
-
1665  AMREX_ALWAYS_ASSERT(real_width >= 0);
-
1666  AMREX_ALWAYS_ASSERT(real_set_width >= 0);
-
1667  AMREX_ALWAYS_ASSERT(real_width >= real_set_width);
-
1668 
-
1669  if (cf_width < 0 || cf_set_width < 0 || cf_width < cf_set_width) {
-
1670  Abort("You must set cf_width >= cf_set_width >= 0");
-
1671  }
-
1672  if (max_level > 0 && cf_set_width > 0) {
-
1673  for (int lev = 1; lev <= max_level; lev++) {
-
1674  if (cf_set_width%ref_ratio[lev-1][0] != 0 ||
-
1675  cf_set_width%ref_ratio[lev-1][1] != 0 ||
-
1676  cf_set_width%ref_ratio[lev-1][2] != 0 ) {
-
1677  Abort("You must set cf_width to be a multiple of ref_ratio");
-
1678  }
-
1679  }
-
1680  }
+
1676 {
+
1677  AMREX_ALWAYS_ASSERT(cfl > 0. || fixed_dt[0] > 0.);
+
1678 
+
1679  // We don't allow use_real_bcs to be true if init_type is not either InitType::Real or InitType::Metgrid
+
1680  AMREX_ALWAYS_ASSERT(!use_real_bcs || ((init_type == InitType::Real) || (init_type == InitType::Metgrid)) );
1681 
-
1682  // If fixed_mri_dt_ratio is set, it must be even
-
1683  if (fixed_mri_dt_ratio > 0 && (fixed_mri_dt_ratio%2 != 0) )
-
1684  {
-
1685  Abort("If you specify fixed_mri_dt_ratio, it must be even");
-
1686  }
-
1687 
-
1688  for (int lev = 0; lev <= max_level; lev++)
-
1689  {
-
1690  // We ignore fixed_fast_dt if not substepping
-
1691  if (solverChoice.substepping_type[lev] == SubsteppingType::None) {
-
1692  fixed_fast_dt[lev] = -1.0;
-
1693  }
-
1694 
-
1695  // If both fixed_dt and fast_dt are specified, their ratio must be an even integer
-
1696  if (fixed_dt[lev] > 0. && fixed_fast_dt[lev] > 0. && fixed_mri_dt_ratio <= 0)
-
1697  {
-
1698  Real eps = 1.e-12;
-
1699  int ratio = static_cast<int>( ( (1.0+eps) * fixed_dt[lev] ) / fixed_fast_dt[lev] );
-
1700  if (fixed_dt[lev] / fixed_fast_dt[lev] != ratio)
-
1701  {
-
1702  Abort("Ratio of fixed_dt to fixed_fast_dt must be an even integer");
-
1703  }
-
1704  }
-
1705 
-
1706  // If all three are specified, they must be consistent
-
1707  if (fixed_dt[lev] > 0. && fixed_fast_dt[lev] > 0. && fixed_mri_dt_ratio > 0)
-
1708  {
-
1709  if (fixed_dt[lev] / fixed_fast_dt[lev] != fixed_mri_dt_ratio)
-
1710  {
-
1711  Abort("Dt is over-specfied");
-
1712  }
-
1713  }
-
1714  } // lev
-
1715 
-
1716  if (solverChoice.coupling_type == CouplingType::TwoWay && cf_width > 0) {
-
1717  Abort("For two-way coupling you must set cf_width = 0");
-
1718  }
-
1719 }
+
1682  AMREX_ALWAYS_ASSERT(real_width >= 0);
+
1683  AMREX_ALWAYS_ASSERT(real_set_width >= 0);
+
1684  AMREX_ALWAYS_ASSERT(real_width >= real_set_width);
+
1685 
+
1686  if (cf_width < 0 || cf_set_width < 0 || cf_width < cf_set_width) {
+
1687  Abort("You must set cf_width >= cf_set_width >= 0");
+
1688  }
+
1689  if (max_level > 0 && cf_set_width > 0) {
+
1690  for (int lev = 1; lev <= max_level; lev++) {
+
1691  if (cf_set_width%ref_ratio[lev-1][0] != 0 ||
+
1692  cf_set_width%ref_ratio[lev-1][1] != 0 ||
+
1693  cf_set_width%ref_ratio[lev-1][2] != 0 ) {
+
1694  Abort("You must set cf_width to be a multiple of ref_ratio");
+
1695  }
+
1696  }
+
1697  }
+
1698 
+
1699  // If fixed_mri_dt_ratio is set, it must be even
+
1700  if (fixed_mri_dt_ratio > 0 && (fixed_mri_dt_ratio%2 != 0) )
+
1701  {
+
1702  Abort("If you specify fixed_mri_dt_ratio, it must be even");
+
1703  }
+
1704 
+
1705  for (int lev = 0; lev <= max_level; lev++)
+
1706  {
+
1707  // We ignore fixed_fast_dt if not substepping
+
1708  if (solverChoice.substepping_type[lev] == SubsteppingType::None) {
+
1709  fixed_fast_dt[lev] = -1.0;
+
1710  }
+
1711 
+
1712  // If both fixed_dt and fast_dt are specified, their ratio must be an even integer
+
1713  if (fixed_dt[lev] > 0. && fixed_fast_dt[lev] > 0. && fixed_mri_dt_ratio <= 0)
+
1714  {
+
1715  Real eps = 1.e-12;
+
1716  int ratio = static_cast<int>( ( (1.0+eps) * fixed_dt[lev] ) / fixed_fast_dt[lev] );
+
1717  if (fixed_dt[lev] / fixed_fast_dt[lev] != ratio)
+
1718  {
+
1719  Abort("Ratio of fixed_dt to fixed_fast_dt must be an even integer");
+
1720  }
+
1721  }
+
1722 
+
1723  // If all three are specified, they must be consistent
+
1724  if (fixed_dt[lev] > 0. && fixed_fast_dt[lev] > 0. && fixed_mri_dt_ratio > 0)
+
1725  {
+
1726  if (fixed_dt[lev] / fixed_fast_dt[lev] != fixed_mri_dt_ratio)
+
1727  {
+
1728  Abort("Dt is over-specfied");
+
1729  }
+
1730  }
+
1731  } // lev
+
1732 
+
1733  if (solverChoice.coupling_type == CouplingType::TwoWay && cf_width > 0) {
+
1734  Abort("For two-way coupling you must set cf_width = 0");
+
1735  }
+
1736 }
int real_width
Definition: ERF.H:1070
int real_set_width
Definition: ERF.H:1071
@@ -13559,283 +13574,283 @@

-
1378 {
-
1379  {
-
1380  ParmParse pp; // Traditionally, max_step and stop_time do not have prefix.
-
1381  pp.query("max_step", max_step);
-
1382  pp.query("stop_time", stop_time);
-
1383 
-
1384  pp.query("start_time", start_time); // This is optional, it defaults to 0
-
1385  }
-
1386 
-
1387  ParmParse pp(pp_prefix);
-
1388  ParmParse pp_amr("amr");
-
1389  {
-
1390  // The type of the file we restart from
-
1391  pp.query("restart_type", restart_type);
-
1392 
-
1393  pp.query("regrid_int", regrid_int);
-
1394  pp.query("check_file", check_file);
-
1395  pp.query("check_type", check_type);
-
1396 
-
1397  // The regression tests use "amr.restart" and "amr.m_check_int" so we allow
-
1398  // for those or "erf.restart" / "erf.m_check_int" with the former taking
-
1399  // precedence if both are specified
-
1400  pp.query("check_int", m_check_int);
-
1401  pp.query("check_per", m_check_per);
-
1402  pp_amr.query("check_int", m_check_int);
-
1403  pp_amr.query("check_per", m_check_per);
-
1404 
-
1405  pp.query("restart", restart_chkfile);
-
1406  pp_amr.query("restart", restart_chkfile);
-
1407 
-
1408  // Verbosity
-
1409  pp.query("v", verbose);
-
1410  pp.query("mg_v", mg_verbose);
-
1411  pp.query("use_fft", use_fft);
-
1412 #ifndef ERF_USE_FFT
-
1413  if (use_fft) {
-
1414  amrex::Abort("You must build with USE_FFT in order to set use_fft = true in your inputs file");
-
1415  }
-
1416 #endif
-
1417 
-
1418  // Frequency of diagnostic output
-
1419  pp.query("sum_interval", sum_interval);
-
1420  pp.query("sum_period" , sum_per);
+
1395 {
+
1396  {
+
1397  ParmParse pp; // Traditionally, max_step and stop_time do not have prefix.
+
1398  pp.query("max_step", max_step);
+
1399  pp.query("stop_time", stop_time);
+
1400 
+
1401  pp.query("start_time", start_time); // This is optional, it defaults to 0
+
1402  }
+
1403 
+
1404  ParmParse pp(pp_prefix);
+
1405  ParmParse pp_amr("amr");
+
1406  {
+
1407  // The type of the file we restart from
+
1408  pp.query("restart_type", restart_type);
+
1409 
+
1410  pp.query("regrid_int", regrid_int);
+
1411  pp.query("check_file", check_file);
+
1412  pp.query("check_type", check_type);
+
1413 
+
1414  // The regression tests use "amr.restart" and "amr.m_check_int" so we allow
+
1415  // for those or "erf.restart" / "erf.m_check_int" with the former taking
+
1416  // precedence if both are specified
+
1417  pp.query("check_int", m_check_int);
+
1418  pp.query("check_per", m_check_per);
+
1419  pp_amr.query("check_int", m_check_int);
+
1420  pp_amr.query("check_per", m_check_per);
1421 
-
1422  pp.query("pert_interval", pert_interval);
-
1423 
-
1424  // Time step controls
-
1425  pp.query("cfl", cfl);
-
1426  pp.query("substepping_cfl", sub_cfl);
-
1427  pp.query("init_shrink", init_shrink);
-
1428  pp.query("change_max", change_max);
-
1429 
-
1430  fixed_dt.resize(max_level+1,-1.);
-
1431  fixed_fast_dt.resize(max_level+1,-1.);
-
1432 
-
1433  pp.query("fixed_dt", fixed_dt[0]);
-
1434  pp.query("fixed_fast_dt", fixed_fast_dt[0]);
-
1435 
-
1436  for (int lev = 1; lev <= max_level; lev++)
-
1437  {
-
1438  fixed_dt[lev] = fixed_dt[lev-1] / static_cast<Real>(MaxRefRatio(lev-1));
-
1439  fixed_fast_dt[lev] = fixed_fast_dt[lev-1] / static_cast<Real>(MaxRefRatio(lev-1));
-
1440  }
-
1441 
-
1442  pp.query("fixed_mri_dt_ratio", fixed_mri_dt_ratio);
-
1443 
-
1444  // How to initialize
-
1445  init_type = InitType::None;
-
1446  pp.query_enum_case_insensitive("init_type",init_type);
-
1447 
-
1448  // Should we use the bcs we've read in from wrfbdy or metgrid files?
-
1449  // We default to yes if we have them, but the user can override that option
-
1450  use_real_bcs = ( (init_type == InitType::Real) || (init_type == InitType::Metgrid) );
-
1451  pp.query("use_real_bcs",use_real_bcs);
+
1422  pp.query("restart", restart_chkfile);
+
1423  pp_amr.query("restart", restart_chkfile);
+
1424 
+
1425  // Verbosity
+
1426  pp.query("v", verbose);
+
1427  pp.query("mg_v", mg_verbose);
+
1428  pp.query("use_fft", use_fft);
+
1429 #ifndef ERF_USE_FFT
+
1430  if (use_fft) {
+
1431  amrex::Abort("You must build with USE_FFT in order to set use_fft = true in your inputs file");
+
1432  }
+
1433 #endif
+
1434 
+
1435  // Frequency of diagnostic output
+
1436  pp.query("sum_interval", sum_interval);
+
1437  pp.query("sum_period" , sum_per);
+
1438 
+
1439  pp.query("pert_interval", pert_interval);
+
1440 
+
1441  // Time step controls
+
1442  pp.query("cfl", cfl);
+
1443  pp.query("substepping_cfl", sub_cfl);
+
1444  pp.query("init_shrink", init_shrink);
+
1445  pp.query("change_max", change_max);
+
1446 
+
1447  fixed_dt.resize(max_level+1,-1.);
+
1448  fixed_fast_dt.resize(max_level+1,-1.);
+
1449 
+
1450  pp.query("fixed_dt", fixed_dt[0]);
+
1451  pp.query("fixed_fast_dt", fixed_fast_dt[0]);
1452 
-
1453  // We use this to keep track of how many boxes we read in from WRF initialization
-
1454  num_files_at_level.resize(max_level+1,0);
-
1455 
-
1456  // We use this to keep track of how many boxes are specified thru the refinement indicators
-
1457  num_boxes_at_level.resize(max_level+1,0);
-
1458  boxes_at_level.resize(max_level+1);
-
1459 
-
1460  // We always have exactly one file at level 0
-
1461  num_boxes_at_level[0] = 1;
-
1462  boxes_at_level[0].resize(1);
-
1463  boxes_at_level[0][0] = geom[0].Domain();
+
1453  for (int lev = 1; lev <= max_level; lev++)
+
1454  {
+
1455  fixed_dt[lev] = fixed_dt[lev-1] / static_cast<Real>(MaxRefRatio(lev-1));
+
1456  fixed_fast_dt[lev] = fixed_fast_dt[lev-1] / static_cast<Real>(MaxRefRatio(lev-1));
+
1457  }
+
1458 
+
1459  pp.query("fixed_mri_dt_ratio", fixed_mri_dt_ratio);
+
1460 
+
1461  // How to initialize
+
1462  init_type = InitType::None;
+
1463  pp.query_enum_case_insensitive("init_type",init_type);
1464 
-
1465 #ifdef ERF_USE_NETCDF
-
1466  nc_init_file.resize(max_level+1);
-
1467 
-
1468  // NetCDF wrfinput initialization files -- possibly multiple files at each of multiple levels
-
1469  // but we always have exactly one file at level 0
-
1470  for (int lev = 0; lev <= max_level; lev++)
-
1471  {
-
1472  const std::string nc_file_names = Concatenate("nc_init_file_",lev,1);
-
1473  if (pp.contains(nc_file_names.c_str()))
-
1474  {
-
1475  int num_files = pp.countval(nc_file_names.c_str());
-
1476  num_files_at_level[lev] = num_files;
-
1477  nc_init_file[lev].resize(num_files);
-
1478  pp.queryarr(nc_file_names.c_str(), nc_init_file[lev],0,num_files);
-
1479  for (int j = 0; j < num_files; j++)
-
1480  Print() << "Reading NC init file names at level " << lev << " and index " << j << " : " << nc_init_file[lev][j] << std::endl;
-
1481  }
-
1482  }
-
1483 
-
1484  // NetCDF wrfbdy lateral boundary file
-
1485  pp.query("nc_bdy_file", nc_bdy_file);
-
1486 #endif
-
1487 
-
1488  // Flag to trigger initialization from input_sounding like WRF's ideal.exe
-
1489  pp.query("init_sounding_ideal", init_sounding_ideal);
-
1490 
-
1491  // Options for vertical interpolation of met_em*.nc data.
-
1492  pp.query("metgrid_debug_quiescent", metgrid_debug_quiescent);
-
1493  pp.query("metgrid_debug_isothermal", metgrid_debug_isothermal);
-
1494  pp.query("metgrid_debug_dry", metgrid_debug_dry);
-
1495  pp.query("metgrid_debug_psfc", metgrid_debug_psfc);
-
1496  pp.query("metgrid_debug_msf", metgrid_debug_msf);
-
1497  pp.query("metgrid_interp_theta", metgrid_interp_theta);
-
1498  pp.query("metgrid_basic_linear", metgrid_basic_linear);
-
1499  pp.query("metgrid_use_below_sfc", metgrid_use_below_sfc);
-
1500  pp.query("metgrid_use_sfc", metgrid_use_sfc);
-
1501  pp.query("metgrid_retain_sfc", metgrid_retain_sfc);
-
1502  pp.query("metgrid_proximity", metgrid_proximity);
-
1503  pp.query("metgrid_order", metgrid_order);
-
1504  pp.query("metgrid_force_sfc_k", metgrid_force_sfc_k);
-
1505 
-
1506  // Set default to FullState for now ... later we will try Perturbation
-
1507  interpolation_type = StateInterpType::FullState;
-
1508  pp.query_enum_case_insensitive("interpolation_type" ,interpolation_type);
-
1509 
-
1510  PlotFileType plotfile_type_temp = PlotFileType::None;
-
1511  pp.query_enum_case_insensitive("plotfile_type" ,plotfile_type_temp);
-
1512  pp.query_enum_case_insensitive("plotfile_type_1",plotfile_type_1);
-
1513  pp.query_enum_case_insensitive("plotfile_type_2",plotfile_type_2);
-
1514  //
-
1515  // This option is for backward consistency -- if only plotfile_type is set,
-
1516  // then it will be used for both 1 and 2 if and only if they are not set
-
1517  //
-
1518  // Default is native amrex if no type is specified
-
1519  //
-
1520  if (plotfile_type_temp == PlotFileType::None) {
-
1521  if (plotfile_type_1 == PlotFileType::None) {
-
1522  plotfile_type_1 = PlotFileType::Amrex;
-
1523  }
-
1524  if (plotfile_type_2 == PlotFileType::None) {
-
1525  plotfile_type_2 = PlotFileType::Amrex;
-
1526  }
-
1527  } else {
-
1528  if (plotfile_type_1 == PlotFileType::None) {
-
1529  plotfile_type_1 = plotfile_type_temp;
-
1530  } else {
-
1531  amrex::Abort("You must set either plotfile_type or plotfile_type_1, not both");
-
1532  }
-
1533  if (plotfile_type_2 == PlotFileType::None) {
-
1534  plotfile_type_2 = plotfile_type_temp;
-
1535  } else {
-
1536  amrex::Abort("You must set either plotfile_type or plotfile_type_2, not both");
-
1537  }
-
1538  }
-
1539 #ifndef ERF_USE_NETCDF
-
1540  if (plotfile_type_1 == PlotFileType::Netcdf ||
-
1541  plotfile_type_2 == PlotFileType::Netcdf) {
-
1542  amrex::Abort("Plotfile type = Netcdf is not allowed without USE_NETCDF = TRUE");
-
1543  }
-
1544 #endif
-
1545 
-
1546  pp.query("plot_file_1", plot_file_1);
-
1547  pp.query("plot_file_2", plot_file_2);
-
1548  pp.query("plot_int_1" , m_plot_int_1);
-
1549  pp.query("plot_int_2" , m_plot_int_2);
-
1550  pp.query("plot_per_1", m_plot_per_1);
-
1551  pp.query("plot_per_2", m_plot_per_2);
-
1552 
-
1553  pp.query("expand_plotvars_to_unif_rr",m_expand_plotvars_to_unif_rr);
-
1554 
-
1555  pp.query("plot_face_vels",m_plot_face_vels);
-
1556 
-
1557  if ( (m_plot_int_1 > 0 && m_plot_per_1 > 0) ||
-
1558  (m_plot_int_2 > 0 && m_plot_per_2 > 0.) ) {
-
1559  Abort("Must choose only one of plot_int or plot_per");
+
1465  // Should we use the bcs we've read in from wrfbdy or metgrid files?
+
1466  // We default to yes if we have them, but the user can override that option
+
1467  use_real_bcs = ( (init_type == InitType::Real) || (init_type == InitType::Metgrid) );
+
1468  pp.query("use_real_bcs",use_real_bcs);
+
1469 
+
1470  // We use this to keep track of how many boxes we read in from WRF initialization
+
1471  num_files_at_level.resize(max_level+1,0);
+
1472 
+
1473  // We use this to keep track of how many boxes are specified thru the refinement indicators
+
1474  num_boxes_at_level.resize(max_level+1,0);
+
1475  boxes_at_level.resize(max_level+1);
+
1476 
+
1477  // We always have exactly one file at level 0
+
1478  num_boxes_at_level[0] = 1;
+
1479  boxes_at_level[0].resize(1);
+
1480  boxes_at_level[0][0] = geom[0].Domain();
+
1481 
+
1482 #ifdef ERF_USE_NETCDF
+
1483  nc_init_file.resize(max_level+1);
+
1484 
+
1485  // NetCDF wrfinput initialization files -- possibly multiple files at each of multiple levels
+
1486  // but we always have exactly one file at level 0
+
1487  for (int lev = 0; lev <= max_level; lev++)
+
1488  {
+
1489  const std::string nc_file_names = Concatenate("nc_init_file_",lev,1);
+
1490  if (pp.contains(nc_file_names.c_str()))
+
1491  {
+
1492  int num_files = pp.countval(nc_file_names.c_str());
+
1493  num_files_at_level[lev] = num_files;
+
1494  nc_init_file[lev].resize(num_files);
+
1495  pp.queryarr(nc_file_names.c_str(), nc_init_file[lev],0,num_files);
+
1496  for (int j = 0; j < num_files; j++)
+
1497  Print() << "Reading NC init file names at level " << lev << " and index " << j << " : " << nc_init_file[lev][j] << std::endl;
+
1498  }
+
1499  }
+
1500 
+
1501  // NetCDF wrfbdy lateral boundary file
+
1502  pp.query("nc_bdy_file", nc_bdy_file);
+
1503 #endif
+
1504 
+
1505  // Flag to trigger initialization from input_sounding like WRF's ideal.exe
+
1506  pp.query("init_sounding_ideal", init_sounding_ideal);
+
1507 
+
1508  // Options for vertical interpolation of met_em*.nc data.
+
1509  pp.query("metgrid_debug_quiescent", metgrid_debug_quiescent);
+
1510  pp.query("metgrid_debug_isothermal", metgrid_debug_isothermal);
+
1511  pp.query("metgrid_debug_dry", metgrid_debug_dry);
+
1512  pp.query("metgrid_debug_psfc", metgrid_debug_psfc);
+
1513  pp.query("metgrid_debug_msf", metgrid_debug_msf);
+
1514  pp.query("metgrid_interp_theta", metgrid_interp_theta);
+
1515  pp.query("metgrid_basic_linear", metgrid_basic_linear);
+
1516  pp.query("metgrid_use_below_sfc", metgrid_use_below_sfc);
+
1517  pp.query("metgrid_use_sfc", metgrid_use_sfc);
+
1518  pp.query("metgrid_retain_sfc", metgrid_retain_sfc);
+
1519  pp.query("metgrid_proximity", metgrid_proximity);
+
1520  pp.query("metgrid_order", metgrid_order);
+
1521  pp.query("metgrid_force_sfc_k", metgrid_force_sfc_k);
+
1522 
+
1523  // Set default to FullState for now ... later we will try Perturbation
+
1524  interpolation_type = StateInterpType::FullState;
+
1525  pp.query_enum_case_insensitive("interpolation_type" ,interpolation_type);
+
1526 
+
1527  PlotFileType plotfile_type_temp = PlotFileType::None;
+
1528  pp.query_enum_case_insensitive("plotfile_type" ,plotfile_type_temp);
+
1529  pp.query_enum_case_insensitive("plotfile_type_1",plotfile_type_1);
+
1530  pp.query_enum_case_insensitive("plotfile_type_2",plotfile_type_2);
+
1531  //
+
1532  // This option is for backward consistency -- if only plotfile_type is set,
+
1533  // then it will be used for both 1 and 2 if and only if they are not set
+
1534  //
+
1535  // Default is native amrex if no type is specified
+
1536  //
+
1537  if (plotfile_type_temp == PlotFileType::None) {
+
1538  if (plotfile_type_1 == PlotFileType::None) {
+
1539  plotfile_type_1 = PlotFileType::Amrex;
+
1540  }
+
1541  if (plotfile_type_2 == PlotFileType::None) {
+
1542  plotfile_type_2 = PlotFileType::Amrex;
+
1543  }
+
1544  } else {
+
1545  if (plotfile_type_1 == PlotFileType::None) {
+
1546  plotfile_type_1 = plotfile_type_temp;
+
1547  } else {
+
1548  amrex::Abort("You must set either plotfile_type or plotfile_type_1, not both");
+
1549  }
+
1550  if (plotfile_type_2 == PlotFileType::None) {
+
1551  plotfile_type_2 = plotfile_type_temp;
+
1552  } else {
+
1553  amrex::Abort("You must set either plotfile_type or plotfile_type_2, not both");
+
1554  }
+
1555  }
+
1556 #ifndef ERF_USE_NETCDF
+
1557  if (plotfile_type_1 == PlotFileType::Netcdf ||
+
1558  plotfile_type_2 == PlotFileType::Netcdf) {
+
1559  amrex::Abort("Plotfile type = Netcdf is not allowed without USE_NETCDF = TRUE");
1560  }
-
1561 
-
1562  pp.query("profile_int", profile_int);
-
1563  pp.query("destag_profiles", destag_profiles);
-
1564 
-
1565  pp.query("plot_lsm", plot_lsm);
-
1566 #ifdef ERF_USE_RRTMGP
-
1567  pp.query("plot_rad", plot_rad);
-
1568 #endif
+
1561 #endif
+
1562 
+
1563  pp.query("plot_file_1", plot_file_1);
+
1564  pp.query("plot_file_2", plot_file_2);
+
1565  pp.query("plot_int_1" , m_plot_int_1);
+
1566  pp.query("plot_int_2" , m_plot_int_2);
+
1567  pp.query("plot_per_1", m_plot_per_1);
+
1568  pp.query("plot_per_2", m_plot_per_2);
1569 
-
1570  pp.query("output_1d_column", output_1d_column);
-
1571  pp.query("column_per", column_per);
-
1572  pp.query("column_interval", column_interval);
-
1573  pp.query("column_loc_x", column_loc_x);
-
1574  pp.query("column_loc_y", column_loc_y);
-
1575  pp.query("column_file_name", column_file_name);
-
1576 
-
1577  // Sampler output frequency
-
1578  pp.query("sampler_per", sampler_per);
-
1579  pp.query("sampler_interval", sampler_interval);
-
1580 
-
1581  // Specify information about outputting planes of data
-
1582  pp.query("output_bndry_planes", output_bndry_planes);
-
1583  pp.query("bndry_output_planes_interval", bndry_output_planes_interval);
-
1584  pp.query("bndry_output_planes_per", bndry_output_planes_per);
-
1585  pp.query("bndry_output_start_time", bndry_output_planes_start_time);
+
1570  pp.query("expand_plotvars_to_unif_rr",m_expand_plotvars_to_unif_rr);
+
1571 
+
1572  pp.query("plot_face_vels",m_plot_face_vels);
+
1573 
+
1574  if ( (m_plot_int_1 > 0 && m_plot_per_1 > 0) ||
+
1575  (m_plot_int_2 > 0 && m_plot_per_2 > 0.) ) {
+
1576  Abort("Must choose only one of plot_int or plot_per");
+
1577  }
+
1578 
+
1579  pp.query("profile_int", profile_int);
+
1580  pp.query("destag_profiles", destag_profiles);
+
1581 
+
1582  pp.query("plot_lsm", plot_lsm);
+
1583 #ifdef ERF_USE_RRTMGP
+
1584  pp.query("plot_rad", plot_rad);
+
1585 #endif
1586 
-
1587  // Specify whether ingest boundary planes of data
-
1588  pp.query("input_bndry_planes", input_bndry_planes);
-
1589 
-
1590  // Query the set and total widths for wrfbdy interior ghost cells
-
1591  pp.query("real_width", real_width);
-
1592  pp.query("real_set_width", real_set_width);
+
1587  pp.query("output_1d_column", output_1d_column);
+
1588  pp.query("column_per", column_per);
+
1589  pp.query("column_interval", column_interval);
+
1590  pp.query("column_loc_x", column_loc_x);
+
1591  pp.query("column_loc_y", column_loc_y);
+
1592  pp.query("column_file_name", column_file_name);
1593 
-
1594  // Query the set and total widths for crse-fine interior ghost cells
-
1595  pp.query("cf_width", cf_width);
-
1596  pp.query("cf_set_width", cf_set_width);
+
1594  // Sampler output frequency
+
1595  pp.query("sampler_per", sampler_per);
+
1596  pp.query("sampler_interval", sampler_interval);
1597 
-
1598  // AmrMesh iterate on grids?
-
1599  bool iterate(true);
-
1600  pp_amr.query("iterate_grids",iterate);
-
1601  if (!iterate) SetIterateToFalse();
-
1602  }
+
1598  // Specify information about outputting planes of data
+
1599  pp.query("output_bndry_planes", output_bndry_planes);
+
1600  pp.query("bndry_output_planes_interval", bndry_output_planes_interval);
+
1601  pp.query("bndry_output_planes_per", bndry_output_planes_per);
+
1602  pp.query("bndry_output_start_time", bndry_output_planes_start_time);
1603 
-
1604 #ifdef ERF_USE_PARTICLES
-
1605  readTracersParams();
-
1606 #endif
-
1607 
-
1608 #ifdef ERF_USE_MULTIBLOCK
- -
1610 #endif
-
1611 
-
1612  solverChoice.init_params(max_level);
-
1613 
-
1614  // Query the canopy model file name
-
1615  std::string forestfile;
-
1616  solverChoice.do_forest_drag = pp.query("forest_file", forestfile);
- -
1618  for (int lev = 0; lev <= max_level; ++lev) {
-
1619  m_forest_drag[lev] = std::make_unique<ForestDrag>(forestfile);
-
1620  }
-
1621  }
-
1622 
-
1623  // No moving terrain with init real (we must do this after init_params
-
1624  // because that is where we set terrain_type
-
1625  if (init_type == InitType::Real && solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
-
1626  Abort("Moving terrain is not supported with init real");
-
1627  }
+
1604  // Specify whether ingest boundary planes of data
+
1605  pp.query("input_bndry_planes", input_bndry_planes);
+
1606 
+
1607  // Query the set and total widths for wrfbdy interior ghost cells
+
1608  pp.query("real_width", real_width);
+
1609  pp.query("real_set_width", real_set_width);
+
1610 
+
1611  // Query the set and total widths for crse-fine interior ghost cells
+
1612  pp.query("cf_width", cf_width);
+
1613  pp.query("cf_set_width", cf_set_width);
+
1614 
+
1615  // AmrMesh iterate on grids?
+
1616  bool iterate(true);
+
1617  pp_amr.query("iterate_grids",iterate);
+
1618  if (!iterate) SetIterateToFalse();
+
1619  }
+
1620 
+
1621 #ifdef ERF_USE_PARTICLES
+
1622  readTracersParams();
+
1623 #endif
+
1624 
+
1625 #ifdef ERF_USE_MULTIBLOCK
+ +
1627 #endif
1628 
-
1629  // What type of land surface model to use
-
1630  // NOTE: Must be checked after init_params
-
1631  if (solverChoice.lsm_type == LandSurfaceType::SLM) {
-
1632  lsm.SetModel<SLM>();
-
1633  Print() << "SLM land surface model!\n";
-
1634  } else if (solverChoice.lsm_type == LandSurfaceType::MM5) {
-
1635  lsm.SetModel<MM5>();
-
1636  Print() << "MM5 land surface model!\n";
-
1637 #ifdef ERF_USE_NOAH
-
1638  } else if (solverChoice.lsm_type == LandSurfaceType::NOAH) {
-
1639  lsm.SetModel<NOAH>();
-
1640  Print() << "NOAH land surface model!\n";
-
1641 #endif
-
1642  } else if (solverChoice.lsm_type == LandSurfaceType::None) {
-
1643  lsm.SetModel<NullSurf>();
-
1644  Print() << "Null land surface model!\n";
-
1645  } else {
-
1646  Abort("Dont know this LandSurfaceType!") ;
-
1647  }
-
1648 
-
1649  if (verbose > 0) {
-
1650  solverChoice.display(max_level);
-
1651  }
-
1652 
- -
1654 }
+
1629  solverChoice.init_params(max_level);
+
1630 
+
1631  // Query the canopy model file name
+
1632  std::string forestfile;
+
1633  solverChoice.do_forest_drag = pp.query("forest_file", forestfile);
+ +
1635  for (int lev = 0; lev <= max_level; ++lev) {
+
1636  m_forest_drag[lev] = std::make_unique<ForestDrag>(forestfile);
+
1637  }
+
1638  }
+
1639 
+
1640  // No moving terrain with init real (we must do this after init_params
+
1641  // because that is where we set terrain_type
+
1642  if (init_type == InitType::Real && solverChoice.terrain_type == TerrainType::MovingFittedMesh) {
+
1643  Abort("Moving terrain is not supported with init real");
+
1644  }
+
1645 
+
1646  // What type of land surface model to use
+
1647  // NOTE: Must be checked after init_params
+
1648  if (solverChoice.lsm_type == LandSurfaceType::SLM) {
+
1649  lsm.SetModel<SLM>();
+
1650  Print() << "SLM land surface model!\n";
+
1651  } else if (solverChoice.lsm_type == LandSurfaceType::MM5) {
+
1652  lsm.SetModel<MM5>();
+
1653  Print() << "MM5 land surface model!\n";
+
1654 #ifdef ERF_USE_NOAH
+
1655  } else if (solverChoice.lsm_type == LandSurfaceType::NOAH) {
+
1656  lsm.SetModel<NOAH>();
+
1657  Print() << "NOAH land surface model!\n";
+
1658 #endif
+
1659  } else if (solverChoice.lsm_type == LandSurfaceType::None) {
+
1660  lsm.SetModel<NullSurf>();
+
1661  Print() << "Null land surface model!\n";
+
1662  } else {
+
1663  Abort("Dont know this LandSurfaceType!") ;
+
1664  }
+
1665 
+
1666  if (verbose > 0) {
+
1667  solverChoice.display(max_level);
+
1668  }
+
1669 
+ +
1671 }
bool metgrid_basic_linear
Definition: ERF.H:1084
amrex::Vector< amrex::Vector< amrex::Box > > boxes_at_level
Definition: ERF.H:719
bool metgrid_debug_msf
Definition: ERF.H:1082
@@ -13857,7 +13872,7 @@

bool metgrid_debug_dry
Definition: ERF.H:1080

bool metgrid_debug_isothermal
Definition: ERF.H:1079
bool metgrid_debug_psfc
Definition: ERF.H:1081
-
void ParameterSanityChecks()
Definition: ERF.cpp:1658
+
void ParameterSanityChecks()
Definition: ERF.cpp:1675
bool m_expand_plotvars_to_unif_rr
Definition: ERF.H:942
amrex::Vector< int > num_boxes_at_level
Definition: ERF.H:717
std::string check_file
Definition: ERF.H:956
@@ -14273,38 +14288,38 @@

-
545 {
-
546  if (lev > 0 && SolverChoice::mesh_type == MeshType::VariableDz)
-
547  {
-
548  //
-
549  // First interpolate from coarser level
-
550  // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab
-
551  // have been pre-filled - this includes ghost cells both inside and outside
-
552  // the domain
-
553  //
-
554  InterpFromCoarseLevel(*temp_zphys_nd, z_phys_nd[lev]->nGrowVect(),
-
555  IntVect(0,0,0), // do not fill ghost cells outside the domain
-
556  *z_phys_nd[lev-1], 0, 0, 1,
-
557  geom[lev-1], geom[lev],
-
558  refRatio(lev-1), &node_bilinear_interp,
- -
560 
-
561  // This recomputes the fine values using the bottom terrain at the fine resolution,
-
562  // and also fills values of z_phys_nd outside the domain
- +
543 {
+
544  if (lev > 0 && SolverChoice::mesh_type == MeshType::VariableDz)
+
545  {
+
546  //
+
547  // First interpolate from coarser level
+
548  // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab
+
549  // have been pre-filled - this includes ghost cells both inside and outside
+
550  // the domain
+
551  //
+
552  InterpFromCoarseLevel(*temp_zphys_nd, z_phys_nd[lev]->nGrowVect(),
+
553  IntVect(0,0,0), // do not fill ghost cells outside the domain
+
554  *z_phys_nd[lev-1], 0, 0, 1,
+
555  geom[lev-1], geom[lev],
+
556  refRatio(lev-1), &node_bilinear_interp,
+ +
558 
+
559  // This recomputes the fine values using the bottom terrain at the fine resolution,
+
560  // and also fills values of z_phys_nd outside the domain
+ +
562 
+
563  std::swap(temp_zphys_nd, z_phys_nd[lev]);
564 
-
565  std::swap(temp_zphys_nd, z_phys_nd[lev]);
+
565  } // use_terrain && lev > 0
566 
-
567  } // use_terrain && lev > 0
-
568 
-
569  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+
567  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
+
568  //
+
569  // This assumes we have already remade the EBGeometry
570  //
-
571  // This assumes we have already remade the EBGeometry
-
572  //
-
573  terrain_blanking[lev]->setVal(1.0);
-
574  MultiFab::Subtract(*terrain_blanking[lev], EBFactory(lev).getVolFrac(), 0, 0, 1, 0);
-
575  }
-
576 }
+
571  terrain_blanking[lev]->setVal(1.0);
+
572  MultiFab::Subtract(*terrain_blanking[lev], EBFactory(lev).getVolFrac(), 0, 0, 1, 0);
+
573  }
+
574 }
Here is the call graph for this function:
@@ -14518,7 +14533,7 @@

476  particleData.Redistribute();

477 #endif
478 }
-
void remake_zphys(int lev, amrex::Real time, std::unique_ptr< amrex::MultiFab > &temp_zphys_nd)
Definition: ERF_MakeNewArrays.cpp:544
+
void remake_zphys(int lev, amrex::Real time, std::unique_ptr< amrex::MultiFab > &temp_zphys_nd)
Definition: ERF_MakeNewArrays.cpp:542

@@ -14536,30 +14551,30 @@

-
1235 {
-
1236  // TODO: This could be deleted since ba/dm are not created yet?
-
1237  for (int lev = 0; lev <= finest_level; ++lev)
-
1238  {
-
1239  auto& lev_new = vars_new[lev];
-
1240  auto& lev_old = vars_old[lev];
-
1241  lev_new[Vars::cons].setVal(0.); lev_old[Vars::cons].setVal(0.);
-
1242  lev_new[Vars::xvel].setVal(0.); lev_old[Vars::xvel].setVal(0.);
-
1243  lev_new[Vars::yvel].setVal(0.); lev_old[Vars::yvel].setVal(0.);
-
1244  lev_new[Vars::zvel].setVal(0.); lev_old[Vars::zvel].setVal(0.);
-
1245  }
-
1246 
-
1247 #ifdef ERF_USE_NETCDF
-
1248  if (restart_type == "netcdf") {
-
1249  ReadNCCheckpointFile();
-
1250  }
-
1251 #endif
-
1252  if (restart_type == "native") {
- -
1254  }
-
1255 
-
1256  // We set this here so that we don't over-write the checkpoint file we just started from
- -
1258 }
+
1252 {
+
1253  // TODO: This could be deleted since ba/dm are not created yet?
+
1254  for (int lev = 0; lev <= finest_level; ++lev)
+
1255  {
+
1256  auto& lev_new = vars_new[lev];
+
1257  auto& lev_old = vars_old[lev];
+
1258  lev_new[Vars::cons].setVal(0.); lev_old[Vars::cons].setVal(0.);
+
1259  lev_new[Vars::xvel].setVal(0.); lev_old[Vars::xvel].setVal(0.);
+
1260  lev_new[Vars::yvel].setVal(0.); lev_old[Vars::yvel].setVal(0.);
+
1261  lev_new[Vars::zvel].setVal(0.); lev_old[Vars::zvel].setVal(0.);
+
1262  }
+
1263 
+
1264 #ifdef ERF_USE_NETCDF
+
1265  if (restart_type == "netcdf") {
+
1266  ReadNCCheckpointFile();
+
1267  }
+
1268 #endif
+
1269  if (restart_type == "native") {
+ +
1271  }
+
1272 
+
1273  // We set this here so that we don't over-write the checkpoint file we just started from
+ +
1275 }
void ReadCheckpointFile()
Definition: ERF_Checkpoint.cpp:273
@@ -16426,94 +16441,94 @@

-
383 {
+
381 {
+
382  // ********************************************************************************************
+
383  // Diffusive terms
384  // ********************************************************************************************
-
385  // Diffusive terms
-
386  // ********************************************************************************************
-
387  bool l_use_terrain = (SolverChoice::terrain_type != TerrainType::None);
-
388  bool l_use_kturb = ( (solverChoice.turbChoice[lev].les_type != LESType::None) ||
-
389  (solverChoice.turbChoice[lev].rans_type != RANSType::None) ||
-
390  (solverChoice.turbChoice[lev].pbl_type != PBLType::None) );
-
391  bool l_use_diff = ( (solverChoice.diffChoice.molec_diff_type != MolecDiffType::None) ||
-
392  l_use_kturb );
-
393  bool l_need_SmnSmn = ( (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) ||
-
394  (solverChoice.turbChoice[lev].rans_type == RANSType::kEqn) );
-
395  bool l_use_moist = ( solverChoice.moisture_type != MoistureType::None );
-
396 
-
397  BoxArray ba12 = convert(ba, IntVect(1,1,0));
-
398  BoxArray ba13 = convert(ba, IntVect(1,0,1));
-
399  BoxArray ba23 = convert(ba, IntVect(0,1,1));
-
400 
-
401  if (l_use_diff) {
-
402  //
-
403  // NOTE: We require ghost cells in the vertical when allowing grids that don't
-
404  // cover the entire vertical extent of the domain at this level
-
405  //
-
406  Tau11_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
-
407  Tau22_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
-
408  Tau33_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
-
409  Tau12_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,1) );
-
410  Tau13_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,1) );
-
411  Tau23_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,1) );
-
412  if (l_use_terrain) {
-
413  Tau21_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,1) );
-
414  Tau31_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,1) );
-
415  Tau32_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,1) );
-
416  } else {
-
417  Tau21_lev[lev] = nullptr;
-
418  Tau31_lev[lev] = nullptr;
-
419  Tau32_lev[lev] = nullptr;
-
420  }
-
421  SFS_hfx1_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(1,0,0)), dm, 1, IntVect(1,1,1) );
-
422  SFS_hfx2_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,1,0)), dm, 1, IntVect(1,1,1) );
-
423  SFS_hfx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
-
424  SFS_diss_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
-
425  SFS_hfx1_lev[lev]->setVal(0.);
-
426  SFS_hfx2_lev[lev]->setVal(0.);
-
427  SFS_hfx3_lev[lev]->setVal(0.);
-
428  SFS_diss_lev[lev]->setVal(0.);
-
429  if (l_use_moist) {
-
430  SFS_q1fx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
-
431  SFS_q2fx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
-
432  SFS_q1fx3_lev[lev]->setVal(0.0);
-
433  SFS_q2fx3_lev[lev]->setVal(0.0);
- -
435  SFS_q1fx1_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(1,0,0)), dm, 1, IntVect(1,1,1) );
-
436  SFS_q1fx2_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,1,0)), dm, 1, IntVect(1,1,1) );
-
437  SFS_q1fx1_lev[lev]->setVal(0.0);
-
438  SFS_q1fx2_lev[lev]->setVal(0.0);
-
439  } else {
-
440  SFS_q1fx1_lev[lev] = nullptr;
-
441  SFS_q1fx2_lev[lev] = nullptr;
-
442  }
-
443  } else {
-
444  SFS_q1fx1_lev[lev] = nullptr;
-
445  SFS_q1fx2_lev[lev] = nullptr;
-
446  SFS_q1fx3_lev[lev] = nullptr;
-
447  SFS_q2fx3_lev[lev] = nullptr;
-
448  }
-
449  } else {
-
450  Tau11_lev[lev] = nullptr; Tau22_lev[lev] = nullptr; Tau33_lev[lev] = nullptr;
-
451  Tau12_lev[lev] = nullptr; Tau21_lev[lev] = nullptr;
-
452  Tau13_lev[lev] = nullptr; Tau31_lev[lev] = nullptr;
-
453  Tau23_lev[lev] = nullptr; Tau32_lev[lev] = nullptr;
-
454  SFS_hfx1_lev[lev] = nullptr; SFS_hfx2_lev[lev] = nullptr; SFS_hfx3_lev[lev] = nullptr;
-
455  SFS_diss_lev[lev] = nullptr;
-
456  }
-
457 
-
458  if (l_use_kturb) {
-
459  eddyDiffs_lev[lev] = std::make_unique<MultiFab>(ba, dm, EddyDiff::NumDiffs, 2);
-
460  eddyDiffs_lev[lev]->setVal(0.0);
-
461  if(l_need_SmnSmn) {
-
462  SmnSmn_lev[lev] = std::make_unique<MultiFab>( ba, dm, 1, 0 );
-
463  } else {
-
464  SmnSmn_lev[lev] = nullptr;
-
465  }
-
466  } else {
-
467  eddyDiffs_lev[lev] = nullptr;
-
468  SmnSmn_lev[lev] = nullptr;
-
469  }
-
470 }
+
385  bool l_use_terrain = (SolverChoice::terrain_type != TerrainType::None);
+
386  bool l_use_kturb = ( (solverChoice.turbChoice[lev].les_type != LESType::None) ||
+
387  (solverChoice.turbChoice[lev].rans_type != RANSType::None) ||
+
388  (solverChoice.turbChoice[lev].pbl_type != PBLType::None) );
+
389  bool l_use_diff = ( (solverChoice.diffChoice.molec_diff_type != MolecDiffType::None) ||
+
390  l_use_kturb );
+
391  bool l_need_SmnSmn = ( (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) ||
+
392  (solverChoice.turbChoice[lev].rans_type == RANSType::kEqn) );
+
393  bool l_use_moist = ( solverChoice.moisture_type != MoistureType::None );
+
394 
+
395  BoxArray ba12 = convert(ba, IntVect(1,1,0));
+
396  BoxArray ba13 = convert(ba, IntVect(1,0,1));
+
397  BoxArray ba23 = convert(ba, IntVect(0,1,1));
+
398 
+
399  if (l_use_diff) {
+
400  //
+
401  // NOTE: We require ghost cells in the vertical when allowing grids that don't
+
402  // cover the entire vertical extent of the domain at this level
+
403  //
+
404  Tau11_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
+
405  Tau22_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
+
406  Tau33_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
+
407  Tau12_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,1) );
+
408  Tau13_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,1) );
+
409  Tau23_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,1) );
+
410  if (l_use_terrain) {
+
411  Tau21_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,1) );
+
412  Tau31_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,1) );
+
413  Tau32_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,1) );
+
414  } else {
+
415  Tau21_lev[lev] = nullptr;
+
416  Tau31_lev[lev] = nullptr;
+
417  Tau32_lev[lev] = nullptr;
+
418  }
+
419  SFS_hfx1_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(1,0,0)), dm, 1, IntVect(1,1,1) );
+
420  SFS_hfx2_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,1,0)), dm, 1, IntVect(1,1,1) );
+
421  SFS_hfx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
+
422  SFS_diss_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
+
423  SFS_hfx1_lev[lev]->setVal(0.);
+
424  SFS_hfx2_lev[lev]->setVal(0.);
+
425  SFS_hfx3_lev[lev]->setVal(0.);
+
426  SFS_diss_lev[lev]->setVal(0.);
+
427  if (l_use_moist) {
+
428  SFS_q1fx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
+
429  SFS_q2fx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
+
430  SFS_q1fx3_lev[lev]->setVal(0.0);
+
431  SFS_q2fx3_lev[lev]->setVal(0.0);
+ +
433  SFS_q1fx1_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(1,0,0)), dm, 1, IntVect(1,1,1) );
+
434  SFS_q1fx2_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,1,0)), dm, 1, IntVect(1,1,1) );
+
435  SFS_q1fx1_lev[lev]->setVal(0.0);
+
436  SFS_q1fx2_lev[lev]->setVal(0.0);
+
437  } else {
+
438  SFS_q1fx1_lev[lev] = nullptr;
+
439  SFS_q1fx2_lev[lev] = nullptr;
+
440  }
+
441  } else {
+
442  SFS_q1fx1_lev[lev] = nullptr;
+
443  SFS_q1fx2_lev[lev] = nullptr;
+
444  SFS_q1fx3_lev[lev] = nullptr;
+
445  SFS_q2fx3_lev[lev] = nullptr;
+
446  }
+
447  } else {
+
448  Tau11_lev[lev] = nullptr; Tau22_lev[lev] = nullptr; Tau33_lev[lev] = nullptr;
+
449  Tau12_lev[lev] = nullptr; Tau21_lev[lev] = nullptr;
+
450  Tau13_lev[lev] = nullptr; Tau31_lev[lev] = nullptr;
+
451  Tau23_lev[lev] = nullptr; Tau32_lev[lev] = nullptr;
+
452  SFS_hfx1_lev[lev] = nullptr; SFS_hfx2_lev[lev] = nullptr; SFS_hfx3_lev[lev] = nullptr;
+
453  SFS_diss_lev[lev] = nullptr;
+
454  }
+
455 
+
456  if (l_use_kturb) {
+
457  eddyDiffs_lev[lev] = std::make_unique<MultiFab>(ba, dm, EddyDiff::NumDiffs, 2);
+
458  eddyDiffs_lev[lev]->setVal(0.0);
+
459  if(l_need_SmnSmn) {
+
460  SmnSmn_lev[lev] = std::make_unique<MultiFab>( ba, dm, 1, 0 );
+
461  } else {
+
462  SmnSmn_lev[lev] = nullptr;
+
463  }
+
464  } else {
+
465  eddyDiffs_lev[lev] = nullptr;
+
466  SmnSmn_lev[lev] = nullptr;
+
467  }
+
468 }
@ NumDiffs
Definition: ERF_IndexDefines.H:162
@@ -16533,14 +16548,14 @@

-
580 {
-
581  if (SolverChoice::mesh_type == MeshType::StretchedDz ||
-
582  SolverChoice::mesh_type == MeshType::VariableDz) {
-
583  make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]);
-
584  make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]);
-
585  make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]);
-
586  }
-
587 }
+
578 {
+
579  if (SolverChoice::mesh_type == MeshType::StretchedDz ||
+
580  SolverChoice::mesh_type == MeshType::VariableDz) {
+
581  make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]);
+
582  make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]);
+
583  make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]);
+
584  }
+
585 }
void make_areas(const Geometry &geom, MultiFab &z_phys_nd, MultiFab &ax, MultiFab &ay, MultiFab &az)
Definition: ERF_TerrainMetrics.cpp:629
void make_J(const Geometry &geom, MultiFab &z_phys_nd, MultiFab &detJ_cc)
Definition: ERF_TerrainMetrics.cpp:591
@@ -18073,47 +18088,47 @@

-
1954 {
-
1955  bool write_now = false;
-
1956 
-
1957  if ( plot_int > 0 && (nstep % plot_int == 0) ) {
-
1958  write_now = true;
-
1959 
-
1960  } else if (plot_per > 0.0) {
-
1961 
-
1962  // Check to see if we've crossed a plot_per interval by comparing
-
1963  // the number of intervals that have elapsed for both the current
-
1964  // time and the time at the beginning of this timestep.
-
1965 
-
1966  const Real eps = std::numeric_limits<Real>::epsilon() * Real(10.0) * std::abs(cur_time);
-
1967 
-
1968  int num_per_old = static_cast<int>(std::floor((cur_time-eps-dt_lev) / plot_per));
-
1969  int num_per_new = static_cast<int>(std::floor((cur_time-eps ) / plot_per));
-
1970 
-
1971  // Before using these, however, we must test for the case where we're
-
1972  // within machine epsilon of the next interval. In that case, increment
-
1973  // the counter, because we have indeed reached the next plot_per interval
-
1974  // at this point.
-
1975 
-
1976  const Real next_plot_time = (num_per_old + 1) * plot_per;
-
1977 
-
1978  if ((num_per_new == num_per_old) && std::abs(cur_time - next_plot_time) <= eps)
-
1979  {
-
1980  num_per_new += 1;
-
1981  }
+
1971 {
+
1972  bool write_now = false;
+
1973 
+
1974  if ( plot_int > 0 && (nstep % plot_int == 0) ) {
+
1975  write_now = true;
+
1976 
+
1977  } else if (plot_per > 0.0) {
+
1978 
+
1979  // Check to see if we've crossed a plot_per interval by comparing
+
1980  // the number of intervals that have elapsed for both the current
+
1981  // time and the time at the beginning of this timestep.
1982 
-
1983  // Similarly, we have to account for the case where the old time is within
-
1984  // machine epsilon of the beginning of this interval, so that we don't double
-
1985  // count that time threshold -- we already plotted at that time on the last timestep.
-
1986 
-
1987  if ((num_per_new != num_per_old) && std::abs((cur_time - dt_lev) - next_plot_time) <= eps)
-
1988  num_per_old += 1;
-
1989 
-
1990  if (num_per_old != num_per_new)
-
1991  write_now = true;
-
1992  }
-
1993  return write_now;
-
1994 }
+
1983  const Real eps = std::numeric_limits<Real>::epsilon() * Real(10.0) * std::abs(cur_time);
+
1984 
+
1985  int num_per_old = static_cast<int>(std::floor((cur_time-eps-dt_lev) / plot_per));
+
1986  int num_per_new = static_cast<int>(std::floor((cur_time-eps ) / plot_per));
+
1987 
+
1988  // Before using these, however, we must test for the case where we're
+
1989  // within machine epsilon of the next interval. In that case, increment
+
1990  // the counter, because we have indeed reached the next plot_per interval
+
1991  // at this point.
+
1992 
+
1993  const Real next_plot_time = (num_per_old + 1) * plot_per;
+
1994 
+
1995  if ((num_per_new == num_per_old) && std::abs(cur_time - next_plot_time) <= eps)
+
1996  {
+
1997  num_per_new += 1;
+
1998  }
+
1999 
+
2000  // Similarly, we have to account for the case where the old time is within
+
2001  // machine epsilon of the beginning of this interval, so that we don't double
+
2002  // count that time threshold -- we already plotted at that time on the last timestep.
+
2003 
+
2004  if ((num_per_new != num_per_old) && std::abs((cur_time - dt_lev) - next_plot_time) <= eps)
+
2005  num_per_old += 1;
+
2006 
+
2007  if (num_per_old != num_per_new)
+
2008  write_now = true;
+
2009  }
+
2010  return write_now;
+
2011 }