ERF Class Reference

#include <ERF.H>

Inheritance diagram for ERF:

Collaboration diagram for ERF:

Public Member Functions
	ERF ()
	~ERF () override
void	ERF_shared ()
	ERF (ERF &&) noexcept=delete
ERF &	operator= (ERF &&other) noexcept=delete
	ERF (const ERF &other)=delete
ERF &	operator= (const ERF &other)=delete
void	Evolve ()
void	ErrorEst (int lev, amrex::TagBoxArray &tags, amrex::Real time, int ngrow) override
void	InitData ()
void	InitData_pre ()
void	InitData_post ()
void	compute_divergence (int lev, amrex::MultiFab &rhs, amrex::Vector< amrex::MultiFab > &mom_mf, amrex::Geometry const &geom_at_lev)
void	project_velocities (int lev, amrex::Real dt, amrex::Vector< amrex::MultiFab > &vars, amrex::MultiFab &p)
void	project_velocities_tb (int lev, amrex::Real dt, amrex::Vector< amrex::MultiFab > &vars, amrex::MultiFab &p)
void	solve_with_gmres (int lev, amrex::Vector< amrex::MultiFab > &rhs, amrex::Vector< amrex::MultiFab > &p, amrex::Vector< amrex::Array< amrex::MultiFab, AMREX_SPACEDIM >> &fluxes)
void	solve_with_mlmg (int lev, amrex::Vector< amrex::MultiFab > &rhs, amrex::Vector< amrex::MultiFab > &p, amrex::Vector< amrex::Array< amrex::MultiFab, AMREX_SPACEDIM >> &fluxes)
amrex::Array< amrex::LinOpBCType, AMREX_SPACEDIM >	get_projection_bc (amrex::Orientation::Side side) const noexcept
bool	projection_has_dirichlet (amrex::Array< amrex::LinOpBCType, AMREX_SPACEDIM > bcs) const
void	init_only (int lev, amrex::Real time)
void	restart ()
bool	writeNow (const amrex::Real cur_time, const amrex::Real dt, const int nstep, const int plot_int, const amrex::Real plot_per)
void	post_timestep (int nstep, amrex::Real time, amrex::Real dt_lev)
void	sum_integrated_quantities (amrex::Real time)
void	write_1D_profiles (amrex::Real time)
void	write_1D_profiles_stag (amrex::Real time)
amrex::Real	cloud_fraction (amrex::Real time)
void	FillBdyCCVels (amrex::Vector< amrex::MultiFab > &mf_cc_vel)
void	sample_points (int lev, amrex::Real time, amrex::IntVect cell, amrex::MultiFab &mf)
void	sample_lines (int lev, amrex::Real time, amrex::IntVect cell, amrex::MultiFab &mf)
void	derive_diag_profiles (amrex::Real time, amrex::Gpu::HostVector< amrex::Real > &h_avg_u, amrex::Gpu::HostVector< amrex::Real > &h_avg_v, amrex::Gpu::HostVector< amrex::Real > &h_avg_w, amrex::Gpu::HostVector< amrex::Real > &h_avg_rho, amrex::Gpu::HostVector< amrex::Real > &h_avg_th, amrex::Gpu::HostVector< amrex::Real > &h_avg_ksgs, amrex::Gpu::HostVector< amrex::Real > &h_avg_Kmv, amrex::Gpu::HostVector< amrex::Real > &h_avg_Khv, amrex::Gpu::HostVector< amrex::Real > &h_avg_qv, amrex::Gpu::HostVector< amrex::Real > &h_avg_qc, amrex::Gpu::HostVector< amrex::Real > &h_avg_qr, amrex::Gpu::HostVector< amrex::Real > &h_avg_wqv, amrex::Gpu::HostVector< amrex::Real > &h_avg_wqc, amrex::Gpu::HostVector< amrex::Real > &h_avg_wqr, amrex::Gpu::HostVector< amrex::Real > &h_avg_qi, amrex::Gpu::HostVector< amrex::Real > &h_avg_qs, amrex::Gpu::HostVector< amrex::Real > &h_avg_qg, amrex::Gpu::HostVector< amrex::Real > &h_avg_uu, amrex::Gpu::HostVector< amrex::Real > &h_avg_uv, amrex::Gpu::HostVector< amrex::Real > &h_avg_uw, amrex::Gpu::HostVector< amrex::Real > &h_avg_vv, amrex::Gpu::HostVector< amrex::Real > &h_avg_vw, amrex::Gpu::HostVector< amrex::Real > &h_avg_ww, amrex::Gpu::HostVector< amrex::Real > &h_avg_uth, amrex::Gpu::HostVector< amrex::Real > &h_avg_vth, amrex::Gpu::HostVector< amrex::Real > &h_avg_wth, amrex::Gpu::HostVector< amrex::Real > &h_avg_thth, amrex::Gpu::HostVector< amrex::Real > &h_avg_ku, amrex::Gpu::HostVector< amrex::Real > &h_avg_kv, amrex::Gpu::HostVector< amrex::Real > &h_avg_kw, amrex::Gpu::HostVector< amrex::Real > &h_avg_p, amrex::Gpu::HostVector< amrex::Real > &h_avg_pu, amrex::Gpu::HostVector< amrex::Real > &h_avg_pv, amrex::Gpu::HostVector< amrex::Real > &h_avg_pw, amrex::Gpu::HostVector< amrex::Real > &h_avg_wthv)
void	derive_diag_profiles_stag (amrex::Real time, amrex::Gpu::HostVector< amrex::Real > &h_avg_u, amrex::Gpu::HostVector< amrex::Real > &h_avg_v, amrex::Gpu::HostVector< amrex::Real > &h_avg_w, amrex::Gpu::HostVector< amrex::Real > &h_avg_rho, amrex::Gpu::HostVector< amrex::Real > &h_avg_th, amrex::Gpu::HostVector< amrex::Real > &h_avg_ksgs, amrex::Gpu::HostVector< amrex::Real > &h_avg_Kmv, amrex::Gpu::HostVector< amrex::Real > &h_avg_Khv, amrex::Gpu::HostVector< amrex::Real > &h_avg_qv, amrex::Gpu::HostVector< amrex::Real > &h_avg_qc, amrex::Gpu::HostVector< amrex::Real > &h_avg_qr, amrex::Gpu::HostVector< amrex::Real > &h_avg_wqv, amrex::Gpu::HostVector< amrex::Real > &h_avg_wqc, amrex::Gpu::HostVector< amrex::Real > &h_avg_wqr, amrex::Gpu::HostVector< amrex::Real > &h_avg_qi, amrex::Gpu::HostVector< amrex::Real > &h_avg_qs, amrex::Gpu::HostVector< amrex::Real > &h_avg_qg, amrex::Gpu::HostVector< amrex::Real > &h_avg_uu, amrex::Gpu::HostVector< amrex::Real > &h_avg_uv, amrex::Gpu::HostVector< amrex::Real > &h_avg_uw, amrex::Gpu::HostVector< amrex::Real > &h_avg_vv, amrex::Gpu::HostVector< amrex::Real > &h_avg_vw, amrex::Gpu::HostVector< amrex::Real > &h_avg_ww, amrex::Gpu::HostVector< amrex::Real > &h_avg_uth, amrex::Gpu::HostVector< amrex::Real > &h_avg_vth, amrex::Gpu::HostVector< amrex::Real > &h_avg_wth, amrex::Gpu::HostVector< amrex::Real > &h_avg_thth, amrex::Gpu::HostVector< amrex::Real > &h_avg_ku, amrex::Gpu::HostVector< amrex::Real > &h_avg_kv, amrex::Gpu::HostVector< amrex::Real > &h_avg_kw, amrex::Gpu::HostVector< amrex::Real > &h_avg_p, amrex::Gpu::HostVector< amrex::Real > &h_avg_pu, amrex::Gpu::HostVector< amrex::Real > &h_avg_pv, amrex::Gpu::HostVector< amrex::Real > &h_avg_pw, amrex::Gpu::HostVector< amrex::Real > &h_avg_wthv)
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)
void	derive_stress_profiles_stag (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)
amrex::Real	volWgtSumMF (int lev, const amrex::MultiFab &mf, int comp, const amrex::MultiFab &mapfac, bool local, bool finemask)
void	MakeNewLevelFromCoarse (int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
void	RemakeLevel (int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
void	ClearLevel (int lev) override
void	MakeNewLevelFromScratch (int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
amrex::Real	estTimeStep (int lev, long &dt_fast_ratio) const
void	advance_dycore (int level, amrex::Vector< amrex::MultiFab > &state_old, amrex::Vector< amrex::MultiFab > &state_new, amrex::MultiFab &xvel_old, amrex::MultiFab &yvel_old, amrex::MultiFab &zvel_old, amrex::MultiFab &xvel_new, amrex::MultiFab &yvel_new, amrex::MultiFab &zvel_new, amrex::MultiFab &source, amrex::MultiFab &xmom_src, amrex::MultiFab &ymom_src, amrex::MultiFab &zmom_src, amrex::Geometry fine_geom, amrex::Real dt, amrex::Real time)
void	advance_microphysics (int lev, amrex::MultiFab &cons_in, const amrex::Real &dt_advance, const int &iteration, const amrex::Real &time)
void	advance_lsm (int lev, amrex::MultiFab &, const amrex::Real &dt_advance)
amrex::MultiFab &	build_fine_mask (int lev)
void	MakeHorizontalAverages ()
void	MakeDiagnosticAverage (amrex::Vector< amrex::Real > &h_havg, amrex::MultiFab &S, int n)
void	derive_upwp (amrex::Vector< amrex::Real > &h_havg)
void	WritePlotFile (int which, PlotFileType plotfile_type, amrex::Vector< std::string > plot_var_names)
void	WriteMultiLevelPlotfileWithTerrain (const std::string &plotfilename, int nlevels, const amrex::Vector< const amrex::MultiFab * > &mf, const amrex::Vector< const amrex::MultiFab * > &mf_nd, const amrex::Vector< std::string > &varnames, const amrex::Vector< amrex::Geometry > &my_geom, amrex::Real time, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &my_ref_ratio, const std::string &versionName="HyperCLaw-V1.1", const std::string &levelPrefix="Level_", const std::string &mfPrefix="Cell", const amrex::Vector< std::string > &extra_dirs=amrex::Vector< std::string >()) const
void	WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile, int nlevels, const amrex::Vector< amrex::BoxArray > &bArray, const amrex::Vector< std::string > &varnames, const amrex::Vector< amrex::Geometry > &my_geom, amrex::Real time, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &my_ref_ratio, const std::string &versionName, const std::string &levelPrefix, const std::string &mfPrefix) const
void	erf_enforce_hse (int lev, amrex::MultiFab &dens, amrex::MultiFab &pres, amrex::MultiFab &pi, amrex::MultiFab &th, std::unique_ptr< amrex::MultiFab > &z_cc)
void	init_from_input_sounding (int lev)
void	input_sponge (int lev)
void	init_from_hse (int lev)
void	init_immersed_body (int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
void	fill_from_bndryregs (const amrex::Vector< amrex::MultiFab * > &mfs, amrex::Real time)
void	AverageDownTo (int crse_lev, int scomp, int ncomp)
void	WriteCheckpointFile () const
void	ReadCheckpointFile ()
void	ReadCheckpointFileMOST ()
void	writeJobInfo (const std::string &dir) const

Static Public Member Functions
static bool	is_it_time_for_action (int nstep, amrex::Real time, amrex::Real dt, int action_interval, amrex::Real action_per)
static void	writeBuildInfo (std::ostream &os)
static void	print_banner (MPI_Comm, std::ostream &)
static void	print_usage (MPI_Comm, std::ostream &)
static void	print_error (MPI_Comm, const std::string &msg)
static void	print_summary (std::ostream &)
static void	print_tpls (std::ostream &)

Public Attributes
std::string	pp_prefix {"erf"}

Private Member Functions
void	ReadParameters ()
void	ParameterSanityChecks ()
void	AverageDown ()
void	update_diffusive_arrays (int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
void	init_zphys (int lev, amrex::Real time)
void	remake_zphys (int lev, std::unique_ptr< amrex::MultiFab > &temp_zphys_nd)
void	update_terrain_arrays (int lev)
void	Construct_ERFFillPatchers (int lev)
void	Define_ERFFillPatchers (int lev)
void	init1DArrays ()
void	init_bcs ()
void	init_custom (int lev)
void	init_uniform (int lev)
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)
void	turbPert_update (const int lev, const amrex::Real dt)
void	turbPert_amplitude (const int lev)
void	initialize_integrator (int lev, amrex::MultiFab &cons_mf, amrex::MultiFab &vel_mf)
void	make_physbcs (int lev)
void	initializeMicrophysics (const int &)
void	FillPatch (int lev, amrex::Real time, const amrex::Vector< amrex::MultiFab * > &mfs_vel, bool cons_only=false)
void	FillPatch (int lev, amrex::Real time, const amrex::Vector< amrex::MultiFab * > &mfs_vel, const amrex::Vector< amrex::MultiFab * > &mfs_mom, const amrex::MultiFab &old_base_state, const amrex::MultiFab &new_base_state, bool fillset=true, bool cons_only=false)
void	FillIntermediatePatch (int lev, amrex::Real time, const amrex::Vector< amrex::MultiFab * > &mfs_vel, const amrex::Vector< amrex::MultiFab * > &mfs_mom, int ng_cons, int ng_vel, bool cons_only, int icomp_cons, int ncomp_cons, bool allow_most_bcs=true)
void	FillCoarsePatch (int lev, amrex::Real time)
void	timeStep (int lev, amrex::Real time, int iteration)
void	Advance (int lev, amrex::Real time, amrex::Real dt_lev, int iteration, int ncycle)
void	initHSE ()
	Initialize HSE. More...
void	initHSE (int lev)
void	initRayleigh ()
	Initialize Rayleigh damping profiles. More...
void	initSponge ()
	Initialize sponge profiles. More...
void	setRayleighRefFromSounding (bool restarting)
	Set Rayleigh mean profiles from input sounding. More...
void	setSpongeRefFromSounding (bool restarting)
	Set sponge mean profiles from input sounding. More...
void	ComputeDt (int step=-1)
std::string	PlotFileName (int lev) const
void	setPlotVariables (const std::string &pp_plot_var_names, amrex::Vector< std::string > &plot_var_names)
void	appendPlotVariables (const std::string &pp_plot_var_names, amrex::Vector< std::string > &plot_var_names)
void	init_Dirichlet_bc_data (const std::string input_file)
void	InitializeFromFile ()
void	InitializeLevelFromData (int lev, const amrex::MultiFab &initial_data)
void	post_update (amrex::MultiFab &state_mf, amrex::Real time, const amrex::Geometry &geom)
void	fill_rhs (amrex::MultiFab &rhs_mf, const amrex::MultiFab &state_mf, amrex::Real time, const amrex::Geometry &geom)
void	init_geo_wind_profile (const std::string input_file, amrex::Vector< amrex::Real > &u_geos, amrex::Gpu::DeviceVector< amrex::Real > &u_geos_d, amrex::Vector< amrex::Real > &v_geos, amrex::Gpu::DeviceVector< amrex::Real > &v_geos_d, const amrex::Geometry &lgeom, const amrex::Vector< amrex::Real > &zlev_stag)
void	refinement_criteria_setup ()
AMREX_FORCE_INLINE amrex::YAFluxRegister *	getAdvFluxReg (int lev)
AMREX_FORCE_INLINE std::ostream &	DataLog (int i)
AMREX_FORCE_INLINE int	NumDataLogs () noexcept
AMREX_FORCE_INLINE std::ostream &	SamplePointLog (int i)
AMREX_FORCE_INLINE int	NumSamplePointLogs () noexcept
AMREX_FORCE_INLINE std::ostream &	SampleLineLog (int i)
AMREX_FORCE_INLINE int	NumSampleLineLogs () noexcept
amrex::IntVect &	SamplePoint (int i)
AMREX_FORCE_INLINE int	NumSamplePoints () noexcept
amrex::IntVect &	SampleLine (int i)
AMREX_FORCE_INLINE int	NumSampleLines () noexcept
void	setRecordDataInfo (int i, const std::string &filename)
void	setRecordSamplePointInfo (int i, int lev, amrex::IntVect &cell, const std::string &filename)
void	setRecordSampleLineInfo (int i, int lev, amrex::IntVect &cell, const std::string &filename)
std::string	DataLogName (int i) const noexcept
	The filename of the ith datalog file. More...
std::string	SamplePointLogName (int i) const noexcept
	The filename of the ith sampleptlog file. More...
std::string	SampleLineLogName (int i) const noexcept
	The filename of the ith samplelinelog file. More...

Static Private Member Functions
static amrex::Vector< std::string >	PlotFileVarNames (amrex::Vector< std::string > plot_var_names)
static void	GotoNextLine (std::istream &is)
static AMREX_FORCE_INLINE int	ComputeGhostCells (const AdvChoice &advChoice, bool use_num_diff)
static amrex::Real	getCPUTime ()

Private Attributes
InputSoundingData	input_sounding_data
InputSpongeData	input_sponge_data
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	xvel_bc_data
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	yvel_bc_data
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	zvel_bc_data
std::unique_ptr< ProblemBase >	prob = nullptr
amrex::Vector< int >	num_boxes_at_level
amrex::Vector< int >	num_files_at_level
amrex::Vector< amrex::Vector< amrex::Box > >	boxes_at_level
amrex::Vector< int >	istep
amrex::Vector< int >	nsubsteps
amrex::Vector< amrex::Real >	t_new
amrex::Vector< amrex::Real >	t_old
amrex::Vector< amrex::Real >	dt
amrex::Vector< long >	dt_mri_ratio
amrex::Vector< amrex::Vector< amrex::MultiFab > >	vars_new
amrex::Vector< amrex::Vector< amrex::MultiFab > >	vars_old
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	vel_t_avg
amrex::Vector< amrex::Real >	t_avg_cnt
amrex::Vector< std::unique_ptr< MRISplitIntegrator< amrex::Vector< amrex::MultiFab > > > >	mri_integrator_mem
amrex::Vector< amrex::MultiFab >	pp_inc
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_cons > >	physbcs_cons
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_u > >	physbcs_u
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_v > >	physbcs_v
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_w > >	physbcs_w
amrex::Vector< std::unique_ptr< ERFPhysBCFunct_base > >	physbcs_base
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Theta_prim
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Qv_prim
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Qr_prim
amrex::Vector< amrex::MultiFab >	rU_old
amrex::Vector< amrex::MultiFab >	rU_new
amrex::Vector< amrex::MultiFab >	rV_old
amrex::Vector< amrex::MultiFab >	rV_new
amrex::Vector< amrex::MultiFab >	rW_old
amrex::Vector< amrex::MultiFab >	rW_new
amrex::Vector< amrex::MultiFab >	zmom_crse_rhs
std::unique_ptr< Microphysics >	micro
amrex::Vector< amrex::Vector< amrex::MultiFab * > >	qmoist
LandSurface	lsm
amrex::Vector< amrex::Vector< amrex::MultiFab * > >	lsm_data
amrex::Vector< amrex::Vector< amrex::MultiFab * > >	lsm_flux
int	cf_width {0}
int	cf_set_width {0}
amrex::Vector< ERFFillPatcher >	FPr_c
amrex::Vector< ERFFillPatcher >	FPr_u
amrex::Vector< ERFFillPatcher >	FPr_v
amrex::Vector< ERFFillPatcher >	FPr_w
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau11_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau22_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau33_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau12_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau21_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau13_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau31_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau23_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Tau32_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	eddyDiffs_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SmnSmn_lev
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > >	sst_lev
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::iMultiFab > > >	lmask_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_hfx1_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_hfx2_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_hfx3_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_diss_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_q1fx1_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_q1fx2_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_q1fx3_lev
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	SFS_q2fx3_lev
amrex::Vector< amrex::Vector< amrex::Real > >	zlevels_stag
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	z_phys_nd
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	z_phys_cc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	detJ_cc
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	ax
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	ay
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	az
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	z_phys_nd_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	detJ_cc_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	ax_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	ay_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	az_src
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	z_phys_nd_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	detJ_cc_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	ax_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	ay_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	az_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	z_t_rk
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	mapfac_m
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	mapfac_u
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	mapfac_v
amrex::Vector< amrex::Vector< amrex::Real > >	stretched_dz_h
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	stretched_dz_d
amrex::Vector< amrex::MultiFab >	base_state
amrex::Vector< amrex::MultiFab >	base_state_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Hwave
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Lwave
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Hwave_onegrid
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	Lwave_onegrid
bool	finished_wave = false
amrex::Vector< amrex::YAFluxRegister * >	advflux_reg
amrex::Vector< amrex::BCRec >	domain_bcs_type
amrex::Gpu::DeviceVector< amrex::BCRec >	domain_bcs_type_d
amrex::Array< std::string, 2 *AMREX_SPACEDIM >	domain_bc_type
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max >	m_bc_extdir_vals
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NBCVAR_max >	m_bc_neumann_vals
amrex::GpuArray< ERF_BC, AMREX_SPACEDIM *2 >	phys_bc_type
amrex::Vector< std::unique_ptr< amrex::iMultiFab > >	xflux_imask
amrex::Vector< std::unique_ptr< amrex::iMultiFab > >	yflux_imask
amrex::Vector< std::unique_ptr< amrex::iMultiFab > >	zflux_imask
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	thin_xforce
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	thin_yforce
amrex::Vector< std::unique_ptr< amrex::MultiFab > >	thin_zforce
int	last_plot_file_step_1
int	last_plot_file_step_2
int	last_check_file_step
int	plot_file_on_restart = 1
int	max_step = std::numeric_limits<int>::max()
amrex::Real	start_time = 0.0
amrex::Real	stop_time = std::numeric_limits<amrex::Real>::max()
std::string	restart_chkfile = ""
amrex::Vector< amrex::Real >	fixed_dt
amrex::Vector< amrex::Real >	fixed_fast_dt
int	regrid_int = -1
std::string	plot_file_1 {"plt_1_"}
std::string	plot_file_2 {"plt_2_"}
bool	m_expand_plotvars_to_unif_rr = false
int	m_plot_int_1 = -1
int	m_plot_int_2 = -1
amrex::Real	m_plot_per_1 = -1.0
amrex::Real	m_plot_per_2 = -1.0
bool	m_plot_face_vels = false
bool	plot_lsm = false
int	profile_int = -1
bool	destag_profiles = true
std::string	check_file {"chk"}
std::string	check_type {"native"}
std::string	restart_type {"native"}
int	m_check_int = -1
amrex::Real	m_check_per = -1.0
amrex::Vector< std::string >	plot_var_names_1
amrex::Vector< std::string >	plot_var_names_2
const amrex::Vector< std::string >	cons_names
const amrex::Vector< std::string >	derived_names
TurbulentPerturbation	turbPert
int	real_width {0}
int	real_set_width {0}
amrex::Vector< amrex::Vector< amrex::Real > >	h_rhotheta_src
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	d_rhotheta_src
amrex::Vector< amrex::Vector< amrex::Real > >	h_rhoqt_src
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	d_rhoqt_src
amrex::Vector< amrex::Vector< amrex::Real > >	h_w_subsid
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	d_w_subsid
amrex::Vector< amrex::Vector< amrex::Real > >	h_u_geos
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	d_u_geos
amrex::Vector< amrex::Vector< amrex::Real > >	h_v_geos
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > >	d_v_geos
amrex::Vector< amrex::Vector< amrex::Vector< amrex::Real > > >	h_rayleigh_ptrs
amrex::Vector< amrex::Vector< amrex::Vector< amrex::Real > > >	h_sponge_ptrs
amrex::Vector< amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > >	d_rayleigh_ptrs
amrex::Vector< amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > >	d_sponge_ptrs
amrex::Vector< amrex::Real >	h_havg_density
amrex::Vector< amrex::Real >	h_havg_temperature
amrex::Vector< amrex::Real >	h_havg_pressure
amrex::Vector< amrex::Real >	h_havg_qv
amrex::Vector< amrex::Real >	h_havg_qc
amrex::Gpu::DeviceVector< amrex::Real >	d_havg_density
amrex::Gpu::DeviceVector< amrex::Real >	d_havg_temperature
amrex::Gpu::DeviceVector< amrex::Real >	d_havg_pressure
amrex::Gpu::DeviceVector< amrex::Real >	d_havg_qv
amrex::Gpu::DeviceVector< amrex::Real >	d_havg_qc
std::unique_ptr< WriteBndryPlanes >	m_w2d = nullptr
std::unique_ptr< ReadBndryPlanes >	m_r2d = nullptr
std::unique_ptr< ABLMost >	m_most = nullptr
amrex::Vector< std::unique_ptr< ForestDrag > >	m_forest
amrex::MultiFab	fine_mask
amrex::Vector< amrex::Real >	dz_min
int	sampler_interval = -1
amrex::Real	sampler_per = -1.0
std::unique_ptr< SampleData >	data_sampler = nullptr
amrex::Vector< std::unique_ptr< std::fstream > >	datalog
amrex::Vector< std::string >	datalogname
amrex::Vector< std::unique_ptr< std::fstream > >	sampleptlog
amrex::Vector< std::string >	sampleptlogname
amrex::Vector< amrex::IntVect >	samplepoint
amrex::Vector< std::unique_ptr< std::fstream > >	samplelinelog
amrex::Vector< std::string >	samplelinelogname
amrex::Vector< amrex::IntVect >	sampleline
amrex::Vector< std::unique_ptr< amrex::FabFactory< amrex::FArrayBox > > >	m_factory

Static Private Attributes
static amrex::Real	cfl = 0.8
static amrex::Real	init_shrink = 1.0
static amrex::Real	change_max = 1.1
static int	fixed_mri_dt_ratio = 0
static SolverChoice	solverChoice
static int	verbose = 0
static int	mg_verbose = 0
static bool	use_fft = false
static int	sum_interval = -1
static int	pert_interval = -1
static amrex::Real	sum_per = -1.0
static PlotFileType	plotfile_type_1 = PlotFileType::None
static PlotFileType	plotfile_type_2 = PlotFileType::None
static InitType	init_type
static StateInterpType	interpolation_type
static std::string	sponge_type
static bool	use_real_bcs
static amrex::Vector< amrex::Vector< std::string > >	nc_init_file = {{""}}
static std::string	nc_bdy_file
static bool	init_sounding_ideal = false
static int	output_1d_column = 0
static int	column_interval = -1
static amrex::Real	column_per = -1.0
static amrex::Real	column_loc_x = 0.0
static amrex::Real	column_loc_y = 0.0
static std::string	column_file_name = "column_data.nc"
static int	output_bndry_planes = 0
static int	bndry_output_planes_interval = -1
static amrex::Real	bndry_output_planes_per = -1.0
static amrex::Real	bndry_output_planes_start_time = 0.0
static int	input_bndry_planes = 0
static int	ng_dens_hse
static int	ng_pres_hse
static amrex::Vector< amrex::AMRErrorTag >	ref_tags
static amrex::Real	startCPUTime = 0.0
static amrex::Real	previousCPUTimeUsed = 0.0

Detailed Description

Main class in ERF code, instantiated from main.cpp

Constructor & Destructor Documentation

ERF() [1/3]

ERF::ERF

(

)

{
    ERF_shared();
}

~ERF()

ERF::~ERF ( )

overridedefault

ERF() [2/3]

ERF::ERF ( ERF &&  )

deletenoexcept

ERF() [3/3]

ERF::ERF ( const ERF &  other )

delete

Member Function Documentation

Advance()

void ERF::Advance ( int  lev, amrex::Real  time, amrex::Real  dt_lev, int  iteration, int  ncycle  )

private

Function that advances the solution at one level for a single time step – this does some preliminaries then calls erf_advance

Parameters

[in]	lev	level of refinement (coarsest level is 0)
[in]	time	start time for time advance
[in]	dt_lev	time step for this time advance

{
    BL_PROFILE("ERF::Advance()");
 
    // We must swap the pointers so the previous step's "new" is now this step's "old"
    std::swap(vars_old[lev], vars_new[lev]);
 
    MultiFab& S_old = vars_old[lev][Vars::cons];
    MultiFab& S_new = vars_new[lev][Vars::cons];
 
    MultiFab& U_old = vars_old[lev][Vars::xvel];
    MultiFab& V_old = vars_old[lev][Vars::yvel];
    MultiFab& W_old = vars_old[lev][Vars::zvel];
 
    MultiFab& U_new = vars_new[lev][Vars::xvel];
    MultiFab& V_new = vars_new[lev][Vars::yvel];
    MultiFab& W_new = vars_new[lev][Vars::zvel];
 
    // We need to set these because otherwise in the first call to erf_advance we may
    //    read uninitialized data on ghost values in setting the bc's on the velocities
    U_new.setVal(1.e34,U_new.nGrowVect());
    V_new.setVal(1.e34,V_new.nGrowVect());
    W_new.setVal(1.e34,W_new.nGrowVect());
 
    //
    // NOTE: the momenta here are not fillpatched (they are only used as scratch space)
    // If lev == 0 we have already FillPatched this in ERF::TimeStep
    //
//  if (lev == 0) {
//      FillPatch(lev, time, {&S_old, &U_old, &V_old, &W_old});
//  } else {
    if (lev > 0) {
        FillPatch(lev, time, {&S_old, &U_old, &V_old, &W_old},
                             {&S_old, &rU_old[lev], &rV_old[lev], &rW_old[lev]},
                             base_state[lev], base_state[lev]);
    }
 
    //
    // So we must convert the fillpatched to momenta, including the ghost values
    //
    VelocityToMomentum(U_old, rU_old[lev].nGrowVect(),
                       V_old, rV_old[lev].nGrowVect(),
                       W_old, rW_old[lev].nGrowVect(),
                       S_old, rU_old[lev], rV_old[lev], rW_old[lev],
                       Geom(lev).Domain(),
                       domain_bcs_type);
 
    // TODO: Can test on multiple levels later
    // Update the inflow perturbation update time and amplitude
    if (lev == 0) {
        if (solverChoice.pert_type == PerturbationType::Source ||
            solverChoice.pert_type == PerturbationType::Direct)
        {
            turbPert.calc_tpi_update(lev, dt_lev, U_old, V_old, S_old);
        }
 
        // If PerturbationType::Direct is selected, directly add the computed perturbation
        // on the conserved field
        if (solverChoice.pert_type == PerturbationType::Direct)
        {
            auto m_ixtype = S_old.boxArray().ixType(); // Conserved term
            for (MFIter mfi(S_old,TileNoZ()); mfi.isValid(); ++mfi) {
                Box bx  = mfi.tilebox();
                const Array4<Real> &cell_data  = S_old.array(mfi);
                const Array4<const Real> &pert_cell = turbPert.pb_cell.array(mfi);
                turbPert.apply_tpi(lev, bx, RhoTheta_comp, m_ixtype, cell_data, pert_cell);
            }
        }
    }
 
    // configure ABLMost params if used MostWall boundary condition
    if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) {
        if (m_most) {
            IntVect ng = Theta_prim[lev]->nGrowVect();
            MultiFab::Copy(  *Theta_prim[lev], S_old, RhoTheta_comp, 0, 1, ng);
            MultiFab::Divide(*Theta_prim[lev], S_old, Rho_comp     , 0, 1, ng);
            if (solverChoice.moisture_type != MoistureType::None) {
                ng = Qv_prim[lev]->nGrowVect();
 
                MultiFab::Copy(  *Qv_prim[lev], S_old, RhoQ1_comp, 0, 1, ng);
                MultiFab::Divide(*Qv_prim[lev], S_old, Rho_comp  , 0, 1, ng);
 
                if (solverChoice.RhoQr_comp > -1) {
                    MultiFab::Copy(  *Qr_prim[lev], S_old, solverChoice.RhoQr_comp, 0, 1, ng);
                    MultiFab::Divide(*Qr_prim[lev], S_old, Rho_comp  , 0, 1, ng);
                } else {
                    Qr_prim[lev]->setVal(0.0);
                }
            }
            // NOTE: std::swap above causes the field ptrs to be out of date.
            //       Reassign the field ptrs for MAC avg computation.
            m_most->update_mac_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
            m_most->update_pblh(lev, vars_old, z_phys_cc[lev].get(),
                                solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
            m_most->update_fluxes(lev, time);
        }
    }
 
#if defined(ERF_USE_WINDFARM)
    if (solverChoice.windfarm_type != WindFarmType::None) {
        advance_windfarm(Geom(lev), dt_lev, S_old,
                         U_old, V_old, W_old, vars_windfarm[lev], Nturb[lev], SMark[lev], time);
    }
 
#endif
 
    const BoxArray&            ba = S_old.boxArray();
    const DistributionMapping& dm = S_old.DistributionMap();
 
    int nvars = S_old.nComp();
 
    // Source array for conserved cell-centered quantities -- this will be filled
    //     in the call to make_sources in ERF_TI_slow_rhs_fun.H
    MultiFab cc_source(ba,dm,nvars,1); cc_source.setVal(0.0);
 
    // Source arrays for momenta -- these will be filled
    //     in the call to make_mom_sources in ERF_TI_slow_rhs_fun.H
    BoxArray ba_x(ba); ba_x.surroundingNodes(0);
    MultiFab xmom_source(ba_x,dm,nvars,1); xmom_source.setVal(0.0);
 
    BoxArray ba_y(ba); ba_y.surroundingNodes(1);
    MultiFab ymom_source(ba_y,dm,nvars,1); ymom_source.setVal(0.0);
 
    BoxArray ba_z(ba); ba_z.surroundingNodes(2);
    MultiFab zmom_source(ba_z,dm,nvars,1); zmom_source.setVal(0.0);
 
    // We don't need to call FillPatch on cons_mf because we have fillpatch'ed S_old above
    MultiFab cons_mf(ba,dm,nvars,S_old.nGrowVect());
    MultiFab::Copy(cons_mf,S_old,0,0,S_old.nComp(),S_old.nGrowVect());
 
    amrex::Vector<MultiFab> state_old;
    amrex::Vector<MultiFab> state_new;
 
    // **************************************************************************************
    // Here we define state_old and state_new which are to be advanced
    // **************************************************************************************
    // Initial solution
    // Note that "old" and "new" here are relative to each RK stage.
    state_old.push_back(MultiFab(cons_mf    , amrex::make_alias, 0, nvars)); // cons
    state_old.push_back(MultiFab(rU_old[lev], amrex::make_alias, 0,     1)); // xmom
    state_old.push_back(MultiFab(rV_old[lev], amrex::make_alias, 0,     1)); // ymom
    state_old.push_back(MultiFab(rW_old[lev], amrex::make_alias, 0,     1)); // zmom
 
    // Final solution
    // state_new at the end of the last RK stage holds the t^{n+1} data
    state_new.push_back(MultiFab(S_new      , amrex::make_alias, 0, nvars)); // cons
    state_new.push_back(MultiFab(rU_new[lev], amrex::make_alias, 0,     1)); // xmom
    state_new.push_back(MultiFab(rV_new[lev], amrex::make_alias, 0,     1)); // ymom
    state_new.push_back(MultiFab(rW_new[lev], amrex::make_alias, 0,     1)); // zmom
 
    // **************************************************************************************
    // Update the dycore
    // **************************************************************************************
    advance_dycore(lev, state_old, state_new,
                   U_old, V_old, W_old,
                   U_new, V_new, W_new,
                   cc_source, xmom_source, ymom_source, zmom_source,
                   Geom(lev), dt_lev, time);
 
    // **************************************************************************************
    // Update the microphysics (moisture)
    // **************************************************************************************
    advance_microphysics(lev, S_new, dt_lev, iteration, time);
 
    // **************************************************************************************
    // Update the land surface model
    // **************************************************************************************
    advance_lsm(lev, S_new, dt_lev);
 
#if defined(ERF_USE_RRTMGP)
    // **************************************************************************************
    // Update the radiation
    // **************************************************************************************
    advance_radiation(lev, S_new, dt_lev);
#endif
 
#ifdef ERF_USE_PARTICLES
    // **************************************************************************************
    // Update the particle positions
    // **************************************************************************************
   evolveTracers( lev, dt_lev, vars_new, z_phys_nd );
#endif
 
    // ***********************************************************************************************
    // Impose domain boundary conditions here so that in FillPatching the fine data we won't
    // need to re-fill these
    // ***********************************************************************************************
    if (lev < finest_level) {
         IntVect ngvect_vels = vars_new[lev][Vars::xvel].nGrowVect();
         (*physbcs_cons[lev])(vars_new[lev][Vars::cons],0,vars_new[lev][Vars::cons].nComp(),
                              vars_new[lev][Vars::cons].nGrowVect(),time,BCVars::cons_bc,true);
            (*physbcs_u[lev])(vars_new[lev][Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc,true);
            (*physbcs_v[lev])(vars_new[lev][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc,true);
            (*physbcs_w[lev])(vars_new[lev][Vars::zvel], vars_new[lev][Vars::xvel], vars_new[lev][Vars::yvel],
                              ngvect_vels,time,BCVars::zvel_bc,true);
    }
 
    // **************************************************************************************
    // Register old and new coarse data if we are at a level less than the finest level
    // **************************************************************************************
    if (lev < finest_level)
    {
        if (cf_width > 0) {
            // We must fill the ghost cells of these so that the parallel copy works correctly
            state_old[IntVars::cons].FillBoundary(geom[lev].periodicity());
            state_new[IntVars::cons].FillBoundary(geom[lev].periodicity());
            FPr_c[lev].RegisterCoarseData({&state_old[IntVars::cons], &state_new[IntVars::cons]},
                                          {time, time + dt_lev});
        }
 
        if (cf_width >= 0) {
            // We must fill the ghost cells of these so that the parallel copy works correctly
            state_old[IntVars::xmom].FillBoundary(geom[lev].periodicity());
            state_new[IntVars::xmom].FillBoundary(geom[lev].periodicity());
            FPr_u[lev].RegisterCoarseData({&state_old[IntVars::xmom], &state_new[IntVars::xmom]},
                                          {time, time + dt_lev});
 
            state_old[IntVars::ymom].FillBoundary(geom[lev].periodicity());
            state_new[IntVars::ymom].FillBoundary(geom[lev].periodicity());
            FPr_v[lev].RegisterCoarseData({&state_old[IntVars::ymom], &state_new[IntVars::ymom]},
                                          {time, time + dt_lev});
 
            state_old[IntVars::zmom].FillBoundary(geom[lev].periodicity());
            state_new[IntVars::zmom].FillBoundary(geom[lev].periodicity());
            FPr_w[lev].RegisterCoarseData({&state_old[IntVars::zmom], &state_new[IntVars::zmom]},
                                          {time, time + dt_lev});
        }
 
            //
            // Now create a MultiFab that holds (S_new - S_old) / dt from the coarse level interpolated
            //     on to the coarse/fine boundary at the fine resolution
            //
            Interpolater* mapper_f = &face_cons_linear_interp;
 
            // PhysBCFunctNoOp null_bc;
            // MultiFab tempx(vars_new[lev+1][Vars::xvel].boxArray(),vars_new[lev+1][Vars::xvel].DistributionMap(),1,0);
            // tempx.setVal(0.0);
            // xmom_crse_rhs[lev+1].setVal(0.0);
            // FPr_u[lev].FillSet(tempx               , time       , null_bc, domain_bcs_type);
            // FPr_u[lev].FillSet(xmom_crse_rhs[lev+1], time+dt_lev, null_bc, domain_bcs_type);
            // MultiFab::Subtract(xmom_crse_rhs[lev+1],tempx,0,0,1,IntVect{0});
            // xmom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);
 
            // MultiFab tempy(vars_new[lev+1][Vars::yvel].boxArray(),vars_new[lev+1][Vars::yvel].DistributionMap(),1,0);
            // tempy.setVal(0.0);
            // ymom_crse_rhs[lev+1].setVal(0.0);
            // FPr_v[lev].FillSet(tempy               , time       , null_bc, domain_bcs_type);
            // FPr_v[lev].FillSet(ymom_crse_rhs[lev+1], time+dt_lev, null_bc, domain_bcs_type);
            // MultiFab::Subtract(ymom_crse_rhs[lev+1],tempy,0,0,1,IntVect{0});
            // ymom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);
 
            MultiFab temp_state(zmom_crse_rhs[lev+1].boxArray(),zmom_crse_rhs[lev+1].DistributionMap(),1,0);
            InterpFromCoarseLevel(temp_state,            IntVect{0}, IntVect{0}, state_old[IntVars::zmom], 0, 0, 1,
                                  geom[lev], geom[lev+1], refRatio(lev), mapper_f, domain_bcs_type, BCVars::zvel_bc);
            InterpFromCoarseLevel(zmom_crse_rhs[lev+1],  IntVect{0}, IntVect{0}, state_new[IntVars::zmom], 0, 0, 1,
                                  geom[lev], geom[lev+1], refRatio(lev), mapper_f, domain_bcs_type, BCVars::zvel_bc);
            MultiFab::Subtract(zmom_crse_rhs[lev+1],temp_state,0,0,1,IntVect{0});
            zmom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);
    }
 
    // ***********************************************************************************************
    // Update the time averaged velocities if they are requested
    // ***********************************************************************************************
    if (solverChoice.time_avg_vel) {
        Time_Avg_Vel_atCC(dt[lev], t_avg_cnt[lev], vel_t_avg[lev].get(), U_new, V_new, W_new);
    }
}

Here is the call graph for this function:

advance_dycore()

void ERF::advance_dycore	(	int	level,
		amrex::Vector< amrex::MultiFab > &	state_old,
		amrex::Vector< amrex::MultiFab > &	state_new,
		amrex::MultiFab &	xvel_old,
		amrex::MultiFab &	yvel_old,
		amrex::MultiFab &	zvel_old,
		amrex::MultiFab &	xvel_new,
		amrex::MultiFab &	yvel_new,
		amrex::MultiFab &	zvel_new,
		amrex::MultiFab &	source,
		amrex::MultiFab &	xmom_src,
		amrex::MultiFab &	ymom_src,
		amrex::MultiFab &	zmom_src,
		amrex::Geometry	fine_geom,
		amrex::Real	dt,
		amrex::Real	time
	)

Function that advances the solution at one level for a single time step – this sets up the multirate time integrator and calls the integrator's advance function

Parameters

[in]	level	level of refinement (coarsest level is 0)
[in]	state_old	old-time conserved variables
[in]	state_new	new-time conserved variables
[in]	xvel_old	old-time x-component of velocity
[in]	yvel_old	old-time y-component of velocity
[in]	zvel_old	old-time z-component of velocity
[in]	xvel_new	new-time x-component of velocity
[in]	yvel_new	new-time y-component of velocity
[in]	zvel_new	new-time z-component of velocity
[in]	cc_src	source term for conserved variables
[in]	xmom_src	source term for x-momenta
[in]	ymom_src	source term for y-momenta
[in]	zmom_src	source term for z-momenta
[in]	fine_geom	container for geometry information at current level
[in]	dt_advance	time step for this time advance
[in]	old_time	old time for this time advance

{
    BL_PROFILE_VAR("erf_advance_dycore()",erf_advance_dycore);
 
    const Box& domain = fine_geom.Domain();
 
    DiffChoice dc    = solverChoice.diffChoice;
    TurbChoice tc    = solverChoice.turbChoice[level];
    SpongeChoice sc  = solverChoice.spongeChoice;
 
    MultiFab r_hse (base_state[level], make_alias, BaseState::r0_comp , 1);
    MultiFab p_hse (base_state[level], make_alias, BaseState::p0_comp , 1);
    MultiFab pi_hse(base_state[level], make_alias, BaseState::pi0_comp, 1);
 
    // These pointers are used in the MRI utility functions
    MultiFab* r0  = &r_hse;
    MultiFab* p0  = &p_hse;
    MultiFab* pi0 = &pi_hse;
 
    Real* dptr_rhotheta_src = solverChoice.custom_rhotheta_forcing ? d_rhotheta_src[level].data() : nullptr;
    Real* dptr_rhoqt_src    = solverChoice.custom_moisture_forcing ? d_rhoqt_src[level].data()    : nullptr;
    Real* dptr_wbar_sub     = solverChoice.custom_w_subsidence     ? d_w_subsid[level].data()     : nullptr;
 
    // Turbulent Perturbation Pointer
    //Real* dptr_rhotheta_src = solverChoice.pert_type ? d_rhotheta_src[level].data() : nullptr;
 
    Vector<Real*> d_rayleigh_ptrs_at_lev;
    d_rayleigh_ptrs_at_lev.resize(Rayleigh::nvars);
    d_rayleigh_ptrs_at_lev[Rayleigh::ubar]     = solverChoice.rayleigh_damp_U   ? d_rayleigh_ptrs[level][Rayleigh::ubar].data() : nullptr;
    d_rayleigh_ptrs_at_lev[Rayleigh::vbar]     = solverChoice.rayleigh_damp_V   ? d_rayleigh_ptrs[level][Rayleigh::vbar].data() : nullptr;
    d_rayleigh_ptrs_at_lev[Rayleigh::wbar]     = solverChoice.rayleigh_damp_W   ? d_rayleigh_ptrs[level][Rayleigh::wbar].data() : nullptr;
    d_rayleigh_ptrs_at_lev[Rayleigh::thetabar] = solverChoice.rayleigh_damp_T   ? d_rayleigh_ptrs[level][Rayleigh::thetabar].data() : nullptr;
 
    Vector<Real*> d_sponge_ptrs_at_lev;
    if(sc.sponge_type=="input_sponge")
    {
        d_sponge_ptrs_at_lev.resize(Sponge::nvars_sponge);
        d_sponge_ptrs_at_lev[Sponge::ubar_sponge]  =  d_sponge_ptrs[level][Sponge::ubar_sponge].data();
        d_sponge_ptrs_at_lev[Sponge::vbar_sponge]  =  d_sponge_ptrs[level][Sponge::vbar_sponge].data();
    }
 
    bool l_use_terrain = (SolverChoice::terrain_type != TerrainType::None);
    bool l_use_diff    = ( (dc.molec_diff_type != MolecDiffType::None) ||
                           (tc.les_type        !=       LESType::None) ||
                           (tc.pbl_type        !=       PBLType::None) );
    bool l_use_kturb   = ( (tc.les_type != LESType::None)   ||
                           (tc.pbl_type != PBLType::None) );
    bool l_use_moisture = ( solverChoice.moisture_type != MoistureType::None );
    bool l_implicit_substepping = ( solverChoice.substepping_type[level] == SubsteppingType::Implicit );
 
    const bool use_most = (m_most != nullptr);
    const bool exp_most = (solverChoice.use_explicit_most);
    amrex::ignore_unused(use_most);
 
    const BoxArray& ba            = state_old[IntVars::cons].boxArray();
    const BoxArray& ba_z          = zvel_old.boxArray();
    const DistributionMapping& dm = state_old[IntVars::cons].DistributionMap();
 
    int num_prim = state_old[IntVars::cons].nComp() - 1;
 
    MultiFab    S_prim  (ba  , dm, num_prim,          state_old[IntVars::cons].nGrowVect());
    MultiFab  pi_stage  (ba  , dm,        1,          state_old[IntVars::cons].nGrowVect());
    MultiFab fast_coeffs(ba_z, dm,        5,          0);
    MultiFab* eddyDiffs = eddyDiffs_lev[level].get();
    MultiFab* SmnSmn    = SmnSmn_lev[level].get();
 
    // **************************************************************************************
    // Compute strain for use in slow RHS, Smagorinsky model, and MOST
    // **************************************************************************************
    {
    BL_PROFILE("erf_advance_strain");
    if (l_use_diff) {
 
        const BCRec* bc_ptr_h = domain_bcs_type.data();
        const GpuArray<Real, AMREX_SPACEDIM> dxInv = fine_geom.InvCellSizeArray();
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
        for ( MFIter mfi(state_new[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
        {
            Box bxcc  = mfi.growntilebox(IntVect(1,1,0));
            Box tbxxy = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
            Box tbxxz = mfi.tilebox(IntVect(1,0,1),IntVect(1,1,0));
            Box tbxyz = mfi.tilebox(IntVect(0,1,1),IntVect(1,1,0));
 
            if (bxcc.smallEnd(2) != domain.smallEnd(2)) {
                 bxcc.growLo(2,1);
                tbxxy.growLo(2,1);
                tbxxz.growLo(2,1);
                tbxyz.growLo(2,1);
            }
 
            if (bxcc.bigEnd(2) != domain.bigEnd(2)) {
                 bxcc.growHi(2,1);
                tbxxy.growHi(2,1);
                tbxxz.growHi(2,1);
                tbxyz.growHi(2,1);
            }
 
            const Array4<const Real> & u = xvel_old.array(mfi);
            const Array4<const Real> & v = yvel_old.array(mfi);
            const Array4<const Real> & w = zvel_old.array(mfi);
 
            Array4<Real> tau11 = Tau11_lev[level].get()->array(mfi);
            Array4<Real> tau22 = Tau22_lev[level].get()->array(mfi);
            Array4<Real> tau33 = Tau33_lev[level].get()->array(mfi);
            Array4<Real> tau12 = Tau12_lev[level].get()->array(mfi);
            Array4<Real> tau13 = Tau13_lev[level].get()->array(mfi);
            Array4<Real> tau23 = Tau23_lev[level].get()->array(mfi);
 
            Array4<Real> tau21  = l_use_terrain ? Tau21_lev[level].get()->array(mfi) : Array4<Real>{};
            Array4<Real> tau31  = l_use_terrain ? Tau31_lev[level].get()->array(mfi) : Array4<Real>{};
            Array4<Real> tau32  = l_use_terrain ? Tau32_lev[level].get()->array(mfi) : Array4<Real>{};
            const Array4<const Real>& z_nd = l_use_terrain ? z_phys_nd[level]->const_array(mfi) : Array4<const Real>{};
 
            const Array4<const Real> mf_m = mapfac_m[level]->array(mfi);
            const Array4<const Real> mf_u = mapfac_u[level]->array(mfi);
            const Array4<const Real> mf_v = mapfac_v[level]->array(mfi);
 
            if (l_use_terrain) {
                ComputeStrain_T(bxcc, tbxxy, tbxxz, tbxyz, domain,
                                u, v, w,
                                tau11, tau22, tau33,
                                tau12, tau13,
                                tau21, tau23,
                                tau31, tau32,
                                z_nd, detJ_cc[level]->const_array(mfi), bc_ptr_h, dxInv,
                                mf_m, mf_u, mf_v);
            } else {
                ComputeStrain_N(bxcc, tbxxy, tbxxz, tbxyz, domain,
                                u, v, w,
                                tau11, tau22, tau33,
                                tau12, tau13, tau23,
                                bc_ptr_h, dxInv,
                                mf_m, mf_u, mf_v);
            }
        } // mfi
    } // l_use_diff
    } // profile
 
#include "ERF_TI_utils.H"
 
    // Additional SFS quantities, calculated once per timestep
    MultiFab* Hfx1 = SFS_hfx1_lev[level].get();
    MultiFab* Hfx2 = SFS_hfx2_lev[level].get();
    MultiFab* Hfx3 = SFS_hfx3_lev[level].get();
    MultiFab* Q1fx1 = SFS_q1fx1_lev[level].get();
    MultiFab* Q1fx2 = SFS_q1fx2_lev[level].get();
    MultiFab* Q1fx3 = SFS_q1fx3_lev[level].get();
    MultiFab* Q2fx3 = SFS_q2fx3_lev[level].get();
    MultiFab* Diss = SFS_diss_lev[level].get();
 
    // *************************************************************************
    // Calculate cell-centered eddy viscosity & diffusivities
    //
    // Notes -- we fill all the data in ghost cells before calling this so
    //    that we can fill the eddy viscosity in the ghost regions and
    //    not have to call a boundary filler on this data itself
    //
    // LES - updates both horizontal and vertical eddy viscosity components
    // PBL - only updates vertical eddy viscosity components so horizontal
    //       components come from the LES model or are left as zero.
    // *************************************************************************
    if (l_use_kturb)
    {
        // NOTE: state_new transfers to state_old for PBL (due to ptr swap in advance)
        const BCRec* bc_ptr_h = domain_bcs_type.data();
        ComputeTurbulentViscosity(xvel_old, yvel_old,
                                  *Tau11_lev[level].get(), *Tau22_lev[level].get(), *Tau33_lev[level].get(),
                                  *Tau12_lev[level].get(), *Tau13_lev[level].get(), *Tau23_lev[level].get(),
                                  state_old[IntVars::cons],
                                  *eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated
                                  fine_geom, *mapfac_u[level], *mapfac_v[level],
                                  z_phys_nd[level], solverChoice,
                                  m_most, exp_most, l_use_moisture, level, bc_ptr_h);
    }
 
    // ***********************************************************************************************
    // Update user-defined source terms -- these are defined once per time step (not per RK stage)
    // ***********************************************************************************************
    if (solverChoice.custom_rhotheta_forcing) {
        prob->update_rhotheta_sources(old_time,
                                      h_rhotheta_src[level], d_rhotheta_src[level],
                                      fine_geom, z_phys_cc[level]);
    }
 
    if (solverChoice.custom_moisture_forcing) {
        prob->update_rhoqt_sources(old_time,
                                   h_rhoqt_src[level], d_rhoqt_src[level],
                                   fine_geom, z_phys_cc[level]);
    }
 
    if (solverChoice.custom_geostrophic_profile) {
        prob->update_geostrophic_profile(old_time,
                                   h_u_geos[level], d_u_geos[level],
                                   h_v_geos[level], d_v_geos[level],
                                   fine_geom, z_phys_cc[level]);
    }
 
    // ***********************************************************************************************
    // Convert old velocity available on faces to old momentum on faces to be used in time integration
    // ***********************************************************************************************
    MultiFab density(state_old[IntVars::cons], make_alias, Rho_comp, 1);
 
    //
    // This is an optimization since we won't need more than one ghost
    // cell of momentum in the integrator if not using NumDiff
    //
    IntVect ngu = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : xvel_old.nGrowVect();
    IntVect ngv = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : yvel_old.nGrowVect();
    IntVect ngw = (!solverChoice.use_NumDiff) ? IntVect(1,1,0) : zvel_old.nGrowVect();
 
    VelocityToMomentum(xvel_old, ngu, yvel_old, ngv, zvel_old, ngw, density,
                       state_old[IntVars::xmom],
                       state_old[IntVars::ymom],
                       state_old[IntVars::zmom],
                       domain, domain_bcs_type);
 
    MultiFab::Copy(xvel_new,xvel_old,0,0,1,xvel_old.nGrowVect());
    MultiFab::Copy(yvel_new,yvel_old,0,0,1,yvel_old.nGrowVect());
    MultiFab::Copy(zvel_new,zvel_old,0,0,1,zvel_old.nGrowVect());
 
    bool fast_only = false;
    bool vel_and_mom_synced = true;
 
    apply_bcs(state_old, old_time,
              state_old[IntVars::cons].nGrow(), state_old[IntVars::xmom].nGrow(),
              fast_only, vel_and_mom_synced);
    cons_to_prim(state_old[IntVars::cons], state_old[IntVars::cons].nGrow());
 
#include "ERF_TI_no_substep_fun.H"
#include "ERF_TI_substep_fun.H"
#include "ERF_TI_slow_rhs_fun.H"
 
    // ***************************************************************************************
    // Setup the integrator and integrate for a single timestep
    // **************************************************************************************
    MRISplitIntegrator<Vector<MultiFab> >& mri_integrator = *mri_integrator_mem[level];
 
    // Define rhs and 'post update' utility function that is called after calculating
    // any state data (e.g. at RK stages or at the end of a timestep)
    mri_integrator.set_slow_rhs_pre(slow_rhs_fun_pre);
    mri_integrator.set_slow_rhs_post(slow_rhs_fun_post);
 
    if (solverChoice.anelastic[level]) {
        mri_integrator.set_slow_rhs_inc(slow_rhs_fun_inc);
    }
 
    mri_integrator.set_fast_rhs(fast_rhs_fun);
    mri_integrator.set_slow_fast_timestep_ratio(fixed_mri_dt_ratio > 0 ? fixed_mri_dt_ratio : dt_mri_ratio[level]);
    mri_integrator.set_no_substep(no_substep_fun);
 
    mri_integrator.advance(state_old, state_new, old_time, dt_advance);
 
    if (verbose) Print() << "Done with advance_dycore at level " << level << std::endl;
}

Here is the call graph for this function:

advance_lsm()

void ERF::advance_lsm	(	int	lev,
		amrex::MultiFab &	,
		const amrex::Real &	dt_advance
	)

{
    if (solverChoice.lsm_type != LandSurfaceType::None) {
        lsm.Advance(lev, dt_advance);
    }
}

advance_microphysics()

void ERF::advance_microphysics	(	int	lev,
		amrex::MultiFab &	cons_in,
		const amrex::Real &	dt_advance,
		const int &	iteration,
		const amrex::Real &	time
	)

{
    if (solverChoice.moisture_type != MoistureType::None) {
        micro->Update_Micro_Vars_Lev(lev, cons);
        micro->Advance(lev, dt_advance, iteration, time, solverChoice, vars_new, z_phys_nd);
        micro->Update_State_Vars_Lev(lev, cons);
    }
}

appendPlotVariables()

void ERF::appendPlotVariables ( const std::string &  pp_plot_var_names, amrex::Vector< std::string > &  plot_var_names  )

private

{
    ParmParse pp(pp_prefix);
 
    Vector<std::string> plot_var_names(0);
    if (pp.contains(pp_plot_var_names.c_str())) {
        std::string nm;
        int nPltVars = pp.countval(pp_plot_var_names.c_str());
        for (int i = 0; i < nPltVars; i++) {
            pp.get(pp_plot_var_names.c_str(), nm, i);
            // Add the named variable to our list of plot variables
            // if it is not already in the list
            if (!containerHasElement(plot_var_names, nm)) {
                plot_var_names.push_back(nm);
            }
        }
    }
 
    Vector<std::string> tmp_plot_names(0);
#ifdef ERF_USE_PARTICLES
    Vector<std::string> particle_mesh_plot_names;
    particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
    for (int i = 0; i < particle_mesh_plot_names.size(); i++) {
        std::string tmp(particle_mesh_plot_names[i]);
        if (containerHasElement(plot_var_names, tmp) ) {
            tmp_plot_names.push_back(tmp);
        }
    }
#endif
 
    for (int i = 0; i < tmp_plot_names.size(); i++) {
        a_plot_var_names.push_back( tmp_plot_names[i] );
    }
 
    // Finally, check to see if we found all the requested variables
    for (const auto& plot_name : plot_var_names) {
        if (!containerHasElement(a_plot_var_names, plot_name)) {
             if (amrex::ParallelDescriptor::IOProcessor()) {
                 Warning("\nWARNING: Requested to plot variable '" + plot_name + "' but it is not available");
             }
        }
    }
}

Here is the call graph for this function:

AverageDown()

void ERF::AverageDown ( )

private

{
    AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay);
 
    int src_comp, num_comp;
    for (int lev = finest_level-1; lev >= 0; --lev)
    {
        // If anelastic we don't average down rho because rho == rho0.
        if (solverChoice.anelastic[lev]) {
            src_comp = 1;
        } else {
            src_comp = 0;
        }
        num_comp = vars_new[0][Vars::cons].nComp() - src_comp;
        AverageDownTo(lev,src_comp,num_comp);
    }
}

AverageDownTo()

void ERF::AverageDownTo	(	int	crse_lev,
		int	scomp,
		int	ncomp
	)

{
    if (solverChoice.anelastic[crse_lev]) {
        AMREX_ALWAYS_ASSERT(scomp == 1);
    } else {
        AMREX_ALWAYS_ASSERT(scomp == 0);
    }
 
    AMREX_ALWAYS_ASSERT(ncomp == vars_new[crse_lev][Vars::cons].nComp() - scomp);
    AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay);
 
    // ******************************************************************************************
    // First do cell-centered quantities
    // The quantity that is conserved is not (rho S), but rather (rho S / m^2) where
    // m is the map scale factor at cell centers
    // Here we pre-divide (rho S) by m^2 before average down
    // ******************************************************************************************
    for (int lev = crse_lev; lev <= crse_lev+1; lev++) {
      for (MFIter mfi(vars_new[lev][Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();
        const Array4<      Real>   cons_arr = vars_new[lev][Vars::cons].array(mfi);
        const Array4<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
        if (SolverChoice::terrain_type != TerrainType::None) {
            const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
            ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
            {
                cons_arr(i,j,k,scomp+n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
            });
        } else {
            ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
            {
                cons_arr(i,j,k,scomp+n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
            });
        }
      } // mfi
    } // lev
 
    int fine_lev = crse_lev+1;
 
    if (interpolation_type == StateInterpType::Perturbational) {
        // Make the fine rho and (rho theta) be perturbational
        MultiFab::Divide(vars_new[fine_lev][Vars::cons],vars_new[fine_lev][Vars::cons],
                         Rho_comp,RhoTheta_comp,1,IntVect{0});
        MultiFab::Subtract(vars_new[fine_lev][Vars::cons],base_state[fine_lev],
                           BaseState::r0_comp,Rho_comp,1,IntVect{0});
        MultiFab::Subtract(vars_new[fine_lev][Vars::cons],base_state[fine_lev],
                           BaseState::th0_comp,RhoTheta_comp,1,IntVect{0});
 
        // Make the crse rho and (rho theta) be perturbational
        MultiFab::Divide(vars_new[crse_lev][Vars::cons],vars_new[crse_lev][Vars::cons],
                         Rho_comp,RhoTheta_comp,1,IntVect{0});
        MultiFab::Subtract(vars_new[crse_lev][Vars::cons],base_state[crse_lev],
                           BaseState::r0_comp,Rho_comp,1,IntVect{0});
        MultiFab::Subtract(vars_new[crse_lev][Vars::cons],base_state[crse_lev],
                           BaseState::th0_comp,RhoTheta_comp,1,IntVect{0});
    }
 
    average_down(vars_new[crse_lev+1][Vars::cons],vars_new[crse_lev  ][Vars::cons],
                 scomp, ncomp, refRatio(crse_lev));
 
    if (interpolation_type == StateInterpType::Perturbational) {
        // Restore the fine data to what it was
        MultiFab::Add(vars_new[fine_lev][Vars::cons],base_state[fine_lev],
                      BaseState::r0_comp,Rho_comp,1,IntVect{0});
        MultiFab::Add(vars_new[fine_lev][Vars::cons],base_state[fine_lev],
                      BaseState::th0_comp,RhoTheta_comp,1,IntVect{0});
        MultiFab::Multiply(vars_new[fine_lev][Vars::cons],vars_new[fine_lev][Vars::cons],
                           Rho_comp,RhoTheta_comp,1,IntVect{0});
 
        // Make the crse data be full values not perturbational
        MultiFab::Add(vars_new[crse_lev][Vars::cons],base_state[crse_lev],
                      BaseState::r0_comp,Rho_comp,1,IntVect{0});
        MultiFab::Add(vars_new[crse_lev][Vars::cons],base_state[crse_lev],
                      BaseState::th0_comp,RhoTheta_comp,1,IntVect{0});
        MultiFab::Multiply(vars_new[crse_lev][Vars::cons],vars_new[crse_lev][Vars::cons],
                           Rho_comp,RhoTheta_comp,1,IntVect{0});
    }
 
    vars_new[crse_lev][Vars::cons].FillBoundary(geom[crse_lev].periodicity());
 
    // Here we multiply (rho S) by m^2 after average down
    for (int lev = crse_lev; lev <= crse_lev+1; lev++) {
      for (MFIter mfi(vars_new[lev][Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();
        const Array4<      Real>   cons_arr = vars_new[lev][Vars::cons].array(mfi);
        const Array4<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
        if (SolverChoice::terrain_type != TerrainType::None) {
            const Array4<const Real>   detJ_arr = detJ_cc[lev]->const_array(mfi);
            ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
            {
                cons_arr(i,j,k,scomp+n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k);
            });
        } else {
            ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
            {
                cons_arr(i,j,k,scomp+n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0));
            });
        }
      } // mfi
    } // lev
 
    // ******************************************************************************************
    // Now average down momenta.
    // Note that vars_new holds velocities not momenta, but we want to do conservative
    //    averaging so we first convert to momentum, then average down, then convert
    //    back to velocities -- only on the valid region
    // ******************************************************************************************
    for (int lev = crse_lev; lev <= crse_lev+1; lev++)
    {
        // FillBoundary for density so we can go back and forth between velocity and momentum
        vars_new[lev][Vars::cons].FillBoundary(geom[lev].periodicity());
 
        VelocityToMomentum(vars_new[lev][Vars::xvel], IntVect(0,0,0),
                           vars_new[lev][Vars::yvel], IntVect(0,0,0),
                           vars_new[lev][Vars::zvel], IntVect(0,0,0),
                           vars_new[lev][Vars::cons],
                             rU_new[lev],
                             rV_new[lev],
                             rW_new[lev],
                           Geom(lev).Domain(),
                           domain_bcs_type);
    }
 
    average_down_faces(rU_new[crse_lev+1], rU_new[crse_lev], refRatio(crse_lev), geom[crse_lev]);
    average_down_faces(rV_new[crse_lev+1], rV_new[crse_lev], refRatio(crse_lev), geom[crse_lev]);
    average_down_faces(rW_new[crse_lev+1], rW_new[crse_lev], refRatio(crse_lev), geom[crse_lev]);
 
    for (int lev = crse_lev; lev <= crse_lev+1; lev++) {
        MomentumToVelocity(vars_new[lev][Vars::xvel],
                           vars_new[lev][Vars::yvel],
                           vars_new[lev][Vars::zvel],
                           vars_new[lev][Vars::cons],
                             rU_new[lev],
                             rV_new[lev],
                             rW_new[lev],
                           Geom(lev).Domain(),
                           domain_bcs_type);
    }
}

Here is the call graph for this function:

build_fine_mask()

MultiFab & ERF::build_fine_mask

(

int

level

)

Helper function for constructing a fine mask, that is, a MultiFab masking coarser data at a lower level by zeroing out covered cells in the fine mask MultiFab we compute.

Parameters

level

Fine level index which masks underlying coarser data

{
    // Mask for zeroing covered cells
    AMREX_ASSERT(level > 0);
 
    const BoxArray& cba = grids[level-1];
    const DistributionMapping& cdm = dmap[level-1];
 
    // TODO -- we should make a vector of these a member of ERF class
    fine_mask.define(cba, cdm, 1, 0, MFInfo());
    fine_mask.setVal(1.0);
 
    BoxArray fba = grids[level];
    iMultiFab ifine_mask = makeFineMask(cba, cdm, fba, ref_ratio[level-1], 1, 0);
 
    const auto  fma =  fine_mask.arrays();
    const auto ifma = ifine_mask.arrays();
    ParallelFor(fine_mask, [=] AMREX_GPU_DEVICE(int bno, int i, int j, int k) noexcept
    {
       fma[bno](i,j,k) = ifma[bno](i,j,k);
    });
 
    return fine_mask;
}

ClearLevel()

void ERF::ClearLevel ( int  lev )

override

{
    for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx) {
        vars_new[lev][var_idx].clear();
        vars_old[lev][var_idx].clear();
    }
 
    base_state[lev].clear();
 
    rU_new[lev].clear();
    rU_old[lev].clear();
    rV_new[lev].clear();
    rV_old[lev].clear();
    rW_new[lev].clear();
    rW_old[lev].clear();
 
    if (lev > 0) {
        zmom_crse_rhs[lev].clear();
    }
 
    if (solverChoice.anelastic[lev] == 1) {
        pp_inc[lev].clear();
    }
 
    // Clears the integrator memory
    mri_integrator_mem[lev].reset();
 
    // Clears the physical boundary condition routines
    physbcs_cons[lev].reset();
    physbcs_u[lev].reset();
    physbcs_v[lev].reset();
    physbcs_w[lev].reset();
    physbcs_base[lev].reset();
 
    // Clears the flux register array
    advflux_reg[lev]->reset();
}

cloud_fraction()

Real ERF::cloud_fraction

(

amrex::Real

time

)

{
    BL_PROFILE("ERF::cloud_fraction()");
 
    int lev = 0;
    // This holds all of qc
    MultiFab qc(vars_new[lev][Vars::cons],make_alias,RhoQ2_comp,1);
 
    int direction = 2; // z-direction
    Box const& domain = geom[lev].Domain();
 
    auto const& qc_arr = qc.const_arrays();
 
    // qc_2d is an BaseFab<int> holding the max value over the column
    auto qc_2d = ReduceToPlane<ReduceOpMax,int>(direction, domain, qc,
         [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> int
         {
             if (qc_arr[box_no](i,j,k) > 0) {
                 return 1;
             } else {
                 return 0;
             }
         });
 
    auto* p = qc_2d.dataPtr();
 
    Long numpts = qc_2d.numPts();
 
    AMREX_ASSERT(numpts < Long(std::numeric_limits<int>::max));
 
#if 1
    if (ParallelDescriptor::UseGpuAwareMpi()) {
        ParallelDescriptor::ReduceIntMax(p,static_cast<int>(numpts));
    } else {
        Gpu::PinnedVector<int> hv(numpts);
        Gpu::copyAsync(Gpu::deviceToHost, p, p+numpts, hv.data());
        Gpu::streamSynchronize();
        ParallelDescriptor::ReduceIntMax(hv.data(),static_cast<int>(numpts));
        Gpu::copyAsync(Gpu::hostToDevice, hv.data(), hv.data()+numpts, p);
    }
 
    // Sum over component 0
    Long num_cloudy = qc_2d.template sum<RunOn::Device>(0);
 
#else
    //
    // We need this if we allow domain decomposition in the vertical
    //    but for now we leave it commented out
    //
    Long num_cloudy = Reduce::Sum<Long>(numpts,
         [=] AMREX_GPU_DEVICE (Long i) -> Long {
             if (p[i] == 1) {
                 return 1;
             } else {
                 return 0;
             }
         });
    ParallelDescriptor::ReduceLongSum(num_cloudy);
#endif
 
    Real num_total = qc_2d.box().d_numPts();
 
    Real cloud_frac = num_cloudy / num_total;
 
    return cloud_frac;
}

compute_divergence()

void ERF::compute_divergence	(	int	lev,
		amrex::MultiFab &	rhs,
		amrex::Vector< amrex::MultiFab > &	mom_mf,
		amrex::Geometry const &	geom_at_lev
	)

Project the single-level velocity field to enforce incompressibility Note that the level may or may not be level 0.

{
    BL_PROFILE("ERF::compute_divergence()");
 
    bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None);
 
    auto dxInv = geom[lev].InvCellSizeArray();
 
    Array<MultiFab const*, AMREX_SPACEDIM> rho0_u_const;
    rho0_u_const[0] = &mom_mf[IntVars::xmom];
    rho0_u_const[1] = &mom_mf[IntVars::ymom];
    rho0_u_const[2] = &mom_mf[IntVars::zmom];
 
    // ****************************************************************************
    // Compute divergence which will form RHS
    // Note that we replace "rho0w" with the contravariant momentum, Omega
    // ****************************************************************************
#ifdef ERF_USE_EB
    bool already_on_centroids = true;
    EB_computeDivergence(rhs, rho0_u_const, geom_at_lev, already_on_centroids);
#else
    if (l_use_terrain && SolverChoice::terrain_is_flat)
    {
        for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            Box bx = mfi.tilebox();
            const Array4<Real const>& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi);
            const Array4<Real const>& rho0v_arr = mom_mf[IntVars::ymom].const_array(mfi);
            const Array4<Real const>& rho0w_arr = mom_mf[IntVars::zmom].const_array(mfi);
            const Array4<Real      >&  rhs_arr = rhs.array(mfi);
 
            Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
                Real dz = stretched_dz_d_ptr[k];
                rhs_arr(i,j,k) =  (rho0u_arr(i+1,j,k) - rho0u_arr(i,j,k)) * dxInv[0]
                                 +(rho0v_arr(i,j+1,k) - rho0v_arr(i,j,k)) * dxInv[1]
                                 +(rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) / dz;
            });
        } // mfi
    }
    else if (l_use_terrain) // terrain is not flat
    {
        //
        // Note we compute the divergence using "rho0w" == Omega
        //
        for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            Box bx = mfi.tilebox();
            const Array4<Real >& rho0u_arr = mom_mf[IntVars::xmom].array(mfi);
            const Array4<Real >& rho0v_arr = mom_mf[IntVars::ymom].array(mfi);
            const Array4<Real >& rho0w_arr = mom_mf[IntVars::zmom].array(mfi);
            const Array4<Real      >&  rhs_arr = rhs.array(mfi);
 
            const Array4<Real const>& ax_arr = ax[lev]->const_array(mfi);
            const Array4<Real const>& ay_arr = ay[lev]->const_array(mfi);
            const Array4<Real const>& az_arr = az[lev]->const_array(mfi);
            const Array4<Real const>& dJ_arr = detJ_cc[lev]->const_array(mfi);
 
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                rhs_arr(i,j,k) =   ((ax_arr(i+1,j,k)*rho0u_arr(i+1,j,k) - ax_arr(i,j,k)*rho0u_arr(i,j,k)) * dxInv[0]
                                   +(ay_arr(i,j+1,k)*rho0v_arr(i,j+1,k) - ay_arr(i,j,k)*rho0v_arr(i,j,k)) * dxInv[1]
                                   +(az_arr(i,j,k+1)*rho0w_arr(i,j,k+1) - az_arr(i,j,k)*rho0w_arr(i,j,k)) * dxInv[2]) / dJ_arr(i,j,k);
            });
        } // mfi
 
    }
    else // no terrain
    {
        computeDivergence(rhs, rho0_u_const, geom_at_lev);
    }
#endif
}

ComputeDt()

void ERF::ComputeDt ( int  step = -1 )

private

Function that calls estTimeStep for each level

{
    Vector<Real> dt_tmp(finest_level+1);
 
    for (int lev = 0; lev <= finest_level; ++lev)
    {
        dt_tmp[lev] = estTimeStep(lev, dt_mri_ratio[lev]);
    }
 
    ParallelDescriptor::ReduceRealMin(&dt_tmp[0], dt_tmp.size());
 
    Real dt_0 = dt_tmp[0];
    int n_factor = 1;
    for (int lev = 0; lev <= finest_level; ++lev) {
        dt_tmp[lev] = amrex::min(dt_tmp[lev], change_max*dt[lev]);
        n_factor *= nsubsteps[lev];
        dt_0 = amrex::min(dt_0, n_factor*dt_tmp[lev]);
        if (step == 0){
            dt_0 *= init_shrink;
            if (verbose && init_shrink != 1.0) {
                Print() << "Timestep 0: shrink initial dt at level " << lev << " by " << init_shrink << std::endl;
            }
        }
    }
 
    // Limit dt's by the value of stop_time.
    const Real eps = 1.e-3*dt_0;
    if (t_new[0] + dt_0 > stop_time - eps) {
        dt_0 = stop_time - t_new[0];
    }
 
    dt[0] = dt_0;
    for (int lev = 1; lev <= finest_level; ++lev) {
        dt[lev] = dt[lev-1] / nsubsteps[lev];
    }
}

ComputeGhostCells()

static AMREX_FORCE_INLINE int ERF::ComputeGhostCells ( const AdvChoice &  advChoice, bool  use_num_diff  )

inlinestaticprivate

    {
        if (use_num_diff)
        {
            return 3;
        } else {
            if (
                 (advChoice.dycore_horiz_adv_type    == AdvType::Centered_6th)
              || (advChoice.dycore_vert_adv_type     == AdvType::Centered_6th)
              || (advChoice.moistscal_horiz_adv_type == AdvType::Centered_6th)
              || (advChoice.moistscal_horiz_adv_type == AdvType::Centered_6th)
              || (advChoice.dryscal_horiz_adv_type   == AdvType::Centered_6th)
              || (advChoice.dryscal_vert_adv_type    == AdvType::Centered_6th) )
            { return 3; }
            else if (
                 (advChoice.dycore_horiz_adv_type    == AdvType::Upwind_5th)
              || (advChoice.dycore_vert_adv_type     == AdvType::Upwind_5th)
              || (advChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th)
              || (advChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th)
              || (advChoice.dryscal_horiz_adv_type   == AdvType::Upwind_5th)
              || (advChoice.dryscal_vert_adv_type    == AdvType::Upwind_5th) )
            { return 3; }
            else if (
                  (advChoice.dryscal_horiz_adv_type   == AdvType::Weno_5)
               || (advChoice.dryscal_vert_adv_type    == AdvType::Weno_5)
               || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_5)
               || (advChoice.moistscal_vert_adv_type  == AdvType::Weno_5)
               || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_5Z)
               || (advChoice.moistscal_vert_adv_type  == AdvType::Weno_5Z)
               || (advChoice.dryscal_horiz_adv_type   == AdvType::Weno_5Z)
               || (advChoice.dryscal_vert_adv_type    == AdvType::Weno_5Z) )
            { return 3; }
            else if (
                  (advChoice.dryscal_horiz_adv_type   == AdvType::Weno_7)
               || (advChoice.dryscal_vert_adv_type    == AdvType::Weno_7)
               || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_7)
               || (advChoice.moistscal_vert_adv_type  == AdvType::Weno_7)
               || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_7Z)
               || (advChoice.moistscal_vert_adv_type  == AdvType::Weno_7Z)
               || (advChoice.dryscal_horiz_adv_type   == AdvType::Weno_7Z)
               || (advChoice.dryscal_vert_adv_type    == AdvType::Weno_7Z) )
            { return 4; }
            else
            { return 2; }
        }
    }

Construct_ERFFillPatchers()

void ERF::Construct_ERFFillPatchers ( int  lev )

private

{
    auto& fine_new = vars_new[lev];
    auto& crse_new = vars_new[lev-1];
    auto& ba_fine  = fine_new[Vars::cons].boxArray();
    auto& ba_crse  = crse_new[Vars::cons].boxArray();
    auto& dm_fine  = fine_new[Vars::cons].DistributionMap();
    auto& dm_crse  = crse_new[Vars::cons].DistributionMap();
 
    int ncomp = vars_new[lev][Vars::cons].nComp();
 
    FPr_c.emplace_back(ba_fine, dm_fine, geom[lev]  ,
                       ba_crse, dm_crse, geom[lev-1],
                       -cf_width, -cf_set_width, ncomp, &cell_cons_interp);
    FPr_u.emplace_back(convert(ba_fine, IntVect(1,0,0)), dm_fine, geom[lev]  ,
                       convert(ba_crse, IntVect(1,0,0)), dm_crse, geom[lev-1],
                       -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
    FPr_v.emplace_back(convert(ba_fine, IntVect(0,1,0)), dm_fine, geom[lev]  ,
                       convert(ba_crse, IntVect(0,1,0)), dm_crse, geom[lev-1],
                       -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
    FPr_w.emplace_back(convert(ba_fine, IntVect(0,0,1)), dm_fine, geom[lev]  ,
                       convert(ba_crse, IntVect(0,0,1)), dm_crse, geom[lev-1],
                       -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
}

DataLog()

AMREX_FORCE_INLINE std::ostream& ERF::DataLog ( int  i )

inlineprivate

    {
        return *datalog[i];
    }

DataLogName()

std::string ERF::DataLogName ( int  i ) const

inlineprivatenoexcept

The filename of the ith datalog file.

{ return datalogname[i]; }

Define_ERFFillPatchers()

void ERF::Define_ERFFillPatchers ( int  lev )

private

{
    auto& fine_new = vars_new[lev];
    auto& crse_new = vars_new[lev-1];
    auto& ba_fine  = fine_new[Vars::cons].boxArray();
    auto& ba_crse  = crse_new[Vars::cons].boxArray();
    auto& dm_fine  = fine_new[Vars::cons].DistributionMap();
    auto& dm_crse  = crse_new[Vars::cons].DistributionMap();
 
    int ncomp = fine_new[Vars::cons].nComp();
 
    FPr_c[lev-1].Define(ba_fine, dm_fine, geom[lev]  ,
                        ba_crse, dm_crse, geom[lev-1],
                        -cf_width, -cf_set_width, ncomp, &cell_cons_interp);
    FPr_u[lev-1].Define(convert(ba_fine, IntVect(1,0,0)), dm_fine, geom[lev]  ,
                        convert(ba_crse, IntVect(1,0,0)), dm_crse, geom[lev-1],
                        -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
    FPr_v[lev-1].Define(convert(ba_fine, IntVect(0,1,0)), dm_fine, geom[lev]  ,
                        convert(ba_crse, IntVect(0,1,0)), dm_crse, geom[lev-1],
                        -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
    FPr_w[lev-1].Define(convert(ba_fine, IntVect(0,0,1)), dm_fine, geom[lev]  ,
                        convert(ba_crse, IntVect(0,0,1)), dm_crse, geom[lev-1],
                        -cf_width, -cf_set_width, 1, &face_cons_linear_interp);
}

derive_diag_profiles()

void ERF::derive_diag_profiles	(	amrex::Real	time,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_u,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_v,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_w,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_rho,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_th,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_ksgs,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_Kmv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_Khv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qc,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qr,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wqv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wqc,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wqr,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qi,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qs,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qg,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uu,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_vv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_vw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_ww,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_vth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_thth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_ku,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_kv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_kw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_p,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_pu,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_pv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_pw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wthv
	)

Computes the profiles for diagnostic quantities.

Parameters

h_avg_u	Profile for x-velocity on Host
h_avg_v	Profile for y-velocity on Host
h_avg_w	Profile for z-velocity on Host
h_avg_rho	Profile for density on Host
h_avg_th	Profile for potential temperature on Host
h_avg_ksgs	Profile for Kinetic Energy on Host
h_avg_uu	Profile for x-velocity squared on Host
h_avg_uv	Profile for x-velocity * y-velocity on Host
h_avg_uw	Profile for x-velocity * z-velocity on Host
h_avg_vv	Profile for y-velocity squared on Host
h_avg_vw	Profile for y-velocity * z-velocity on Host
h_avg_ww	Profile for z-velocity squared on Host
h_avg_uth	Profile for x-velocity * potential temperature on Host
h_avg_uiuiu	Profile for u_i*u_i*u triple product on Host
h_avg_uiuiv	Profile for u_i*u_i*v triple product on Host
h_avg_uiuiw	Profile for u_i*u_i*w triple product on Host
h_avg_p	Profile for pressure perturbation on Host
h_avg_pu	Profile for pressure perturbation * x-velocity on Host
h_avg_pv	Profile for pressure perturbation * y-velocity on Host
h_avg_pw	Profile for pressure perturbation * z-velocity on Host

{
    // We assume that this is always called at level 0
    int lev = 0;
 
    bool l_use_kturb = ((solverChoice.turbChoice[lev].les_type != LESType::None) ||
                        (solverChoice.turbChoice[lev].pbl_type != PBLType::None));
    bool l_use_KE   = ( (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) ||
                        (solverChoice.turbChoice[lev].pbl_type == PBLType::MYNN25) );
 
    // This will hold rho, theta, ksgs, Kmh, Kmv, uu, uv, uw, vv, vw, ww, uth, vth, wth,
    //                  0      1     2    3    4   5   6   7   8   9  10   11   12   13
    //                thth, uiuiu, uiuiv, uiuiw, p, pu, pv, pw, qv, qc, qr, wqv, wqc, wqr,
    //                  14     15     16     17 18  19  20  21  22  23  24   25   26   27
    //                  qi, qs, qg, wthv
    //                  28  29  30    31
    MultiFab mf_out(grids[lev], dmap[lev], 32, 0);
 
    MultiFab mf_vels(grids[lev], dmap[lev], AMREX_SPACEDIM, 0);
 
    MultiFab  u_cc(mf_vels, make_alias, 0, 1); // u at cell centers
    MultiFab  v_cc(mf_vels, make_alias, 1, 1); // v at cell centers
    MultiFab  w_cc(mf_vels, make_alias, 2, 1); // w at cell centers
 
    average_face_to_cellcenter(mf_vels,0,
        Array<const MultiFab*,3>{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]});
 
    int zdir = 2;
    auto domain = geom[0].Domain();
 
    // Sum in the horizontal plane
    h_avg_u  = sumToLine(mf_vels ,0,1,domain,zdir);
    h_avg_v  = sumToLine(mf_vels ,1,1,domain,zdir);
    h_avg_w  = sumToLine(mf_vels ,2,1,domain,zdir);
 
    int hu_size =  h_avg_u.size();
 
    // Divide by the total number of cells we are averaging over
    Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
    for (int k = 0; k < hu_size; ++k) {
        h_avg_u[k] /= area_z; h_avg_v[k] /= area_z; h_avg_w[k] /= area_z;
    }
 
    Gpu::DeviceVector<Real> d_avg_u(hu_size, Real(0.0));
    Gpu::DeviceVector<Real> d_avg_v(hu_size, Real(0.0));
    Gpu::DeviceVector<Real> d_avg_w(hu_size, Real(0.0));
 
#if 0
    auto* avg_u_ptr = d_avg_u.data();
    auto* avg_v_ptr = d_avg_v.data();
    auto* avg_w_ptr = d_avg_w.data();
#endif
 
    Gpu::copy(Gpu::hostToDevice, h_avg_u.begin(), h_avg_u.end(), d_avg_u.begin());
    Gpu::copy(Gpu::hostToDevice, h_avg_v.begin(), h_avg_v.end(), d_avg_v.begin());
    Gpu::copy(Gpu::hostToDevice, h_avg_w.begin(), h_avg_w.end(), d_avg_w.begin());
 
    int nvars = vars_new[lev][Vars::cons].nComp();
    MultiFab mf_cons(vars_new[lev][Vars::cons], make_alias, 0, nvars);
 
    MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1);
 
    bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
 
    for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        const Array4<Real>& fab_arr = mf_out.array(mfi);
        const Array4<Real>& u_cc_arr =  u_cc.array(mfi);
        const Array4<Real>& v_cc_arr =  v_cc.array(mfi);
        const Array4<Real>& w_cc_arr =  w_cc.array(mfi);
        const Array4<Real>& cons_arr = mf_cons.array(mfi);
        const Array4<Real>&   p0_arr = p_hse.array(mfi);
        const Array4<const Real>& eta_arr = (l_use_kturb) ? eddyDiffs_lev[lev]->const_array(mfi) :
                                                            Array4<const Real>{};
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            Real theta = cons_arr(i,j,k,RhoTheta_comp) / cons_arr(i,j,k,Rho_comp);
            fab_arr(i, j, k, 0) = cons_arr(i,j,k,Rho_comp);
            fab_arr(i, j, k, 1) = theta;
            Real ksgs = 0.0;
            if (l_use_KE) {
                ksgs = cons_arr(i,j,k,RhoKE_comp) / cons_arr(i,j,k,Rho_comp);
            }
            fab_arr(i, j, k, 2) = ksgs;
#if 1
            if (l_use_kturb) {
                fab_arr(i, j, k, 3) = eta_arr(i,j,k,EddyDiff::Mom_v); // Kmv
                fab_arr(i, j, k, 4) = eta_arr(i,j,k,EddyDiff::Theta_v); // Khv
            } else {
                fab_arr(i, j, k, 3) = 0.0;
                fab_arr(i, j, k, 4) = 0.0;
            }
#else
            // Here we hijack the "Kturb" variable name to print out the resolved kinetic energy
            Real upert = u_cc_arr(i,j,k) - avg_u_ptr[k];
            Real vpert = v_cc_arr(i,j,k) - avg_v_ptr[k];
            Real wpert = w_cc_arr(i,j,k) - avg_w_ptr[k];
            fab_arr(i, j, k, 3) = 0.5 * (upert*upert + vpert*vpert + wpert*wpert);
#endif
            fab_arr(i, j, k, 5) = u_cc_arr(i,j,k) * u_cc_arr(i,j,k);   // u*u
            fab_arr(i, j, k, 6) = u_cc_arr(i,j,k) * v_cc_arr(i,j,k);   // u*v
            fab_arr(i, j, k, 7) = u_cc_arr(i,j,k) * w_cc_arr(i,j,k);   // u*w
            fab_arr(i, j, k, 8) = v_cc_arr(i,j,k) * v_cc_arr(i,j,k);   // v*v
            fab_arr(i, j, k, 9) = v_cc_arr(i,j,k) * w_cc_arr(i,j,k);   // v*w
            fab_arr(i, j, k,10) = w_cc_arr(i,j,k) * w_cc_arr(i,j,k);   // w*w
            fab_arr(i, j, k,11) = u_cc_arr(i,j,k) * theta; // u*th
            fab_arr(i, j, k,12) = v_cc_arr(i,j,k) * theta; // v*th
            fab_arr(i, j, k,13) = w_cc_arr(i,j,k) * theta; // w*th
            fab_arr(i, j, k,14) = theta * theta;           // th*th
 
            // if the number of fields is changed above, then be sure to update
            // the following def!
            Real uiui = fab_arr(i,j,k,5) + fab_arr(i,j,k,8) + fab_arr(i,j,k,10);
            fab_arr(i, j, k,15) = uiui * u_cc_arr(i,j,k);   // (ui*ui)*u
            fab_arr(i, j, k,16) = uiui * v_cc_arr(i,j,k);   // (ui*ui)*v
            fab_arr(i, j, k,17) = uiui * w_cc_arr(i,j,k);   // (ui*ui)*w
 
            if (!use_moisture) {
                Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp));
                p -= p0_arr(i,j,k);
                fab_arr(i, j, k,18) = p;                       // p
                fab_arr(i, j, k,19) = p * u_cc_arr(i,j,k);     // p*u
                fab_arr(i, j, k,20) = p * v_cc_arr(i,j,k);     // p*v
                fab_arr(i, j, k,21) = p * w_cc_arr(i,j,k);     // p*w
                fab_arr(i, j, k,22) = 0.;  // qv
                fab_arr(i, j, k,23) = 0.;  // qc
                fab_arr(i, j, k,24) = 0.;  // qr
                fab_arr(i, j, k,25) = 0.;  // w*qv
                fab_arr(i, j, k,26) = 0.;  // w*qc
                fab_arr(i, j, k,27) = 0.;  // w*qr
                fab_arr(i, j, k,28) = 0.;  // qi
                fab_arr(i, j, k,29) = 0.;  // qs
                fab_arr(i, j, k,30) = 0.;  // qg
                fab_arr(i, j, k,31) = 0.;  // w*thv
            }
        });
    } // mfi
 
    if (use_moisture)
    {
        int n_qstate = micro->Get_Qstate_Size();
 
        for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            const Box& bx = mfi.tilebox();
            const Array4<Real>& fab_arr  = mf_out.array(mfi);
            const Array4<Real>& cons_arr = mf_cons.array(mfi);
            const Array4<Real>& u_cc_arr =  u_cc.array(mfi);
            const Array4<Real>& v_cc_arr =  v_cc.array(mfi);
            const Array4<Real>& w_cc_arr =  w_cc.array(mfi);
            const Array4<Real>&   p0_arr = p_hse.array(mfi);
            const Array4<Real>&   qv_arr = qmoist[0][0]->array(mfi); // TODO: Is this written only on lev 0?
 
            int rhoqr_comp = solverChoice.RhoQr_comp;
 
            ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));
 
                p -= p0_arr(i,j,k);
                fab_arr(i, j, k,18) = p;                       // p
                fab_arr(i, j, k,19) = p * u_cc_arr(i,j,k);     // p*u
                fab_arr(i, j, k,20) = p * v_cc_arr(i,j,k);     // p*v
                fab_arr(i, j, k,21) = p * w_cc_arr(i,j,k);     // p*w
                fab_arr(i, j, k,22) = cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp);  // qv
                fab_arr(i, j, k,23) = cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp);  // qc
                if (rhoqr_comp > -1) {
                    fab_arr(i, j, k,24) = cons_arr(i,j,k,rhoqr_comp) / cons_arr(i,j,k,Rho_comp);  // qr
                } else {
                    fab_arr(i, j, k,24) = Real(0.0);
                }
                fab_arr(i, j, k,25) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp);  // w*qv
                fab_arr(i, j, k,26) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp);  // w*qc
                if (rhoqr_comp > -1) {
                    fab_arr(i, j, k,27) = w_cc_arr(i,j,k) * cons_arr(i,j,k,rhoqr_comp) / cons_arr(i,j,k,Rho_comp);  // w*qr
                } else {
                    fab_arr(i, j, k,27) = Real(0.0);
                }
                if (n_qstate > 3) {
                    fab_arr(i, j, k,28) = cons_arr(i,j,k,RhoQ3_comp) / cons_arr(i,j,k,Rho_comp);  // qi
                    fab_arr(i, j, k,29) = cons_arr(i,j,k,RhoQ5_comp) / cons_arr(i,j,k,Rho_comp);  // qs
                    fab_arr(i, j, k,30) = cons_arr(i,j,k,RhoQ6_comp) / cons_arr(i,j,k,Rho_comp);  // qg
                } else {
                    fab_arr(i, j, k,28) = 0.0;  // qi
                    fab_arr(i, j, k,29) = 0.0;  // qs
                    fab_arr(i, j, k,30) = 0.0;  // qg
                }
                Real ql = fab_arr(i, j, k,22) + fab_arr(i, j, k,23);
                Real theta = cons_arr(i,j,k,RhoTheta_comp) / cons_arr(i,j,k,Rho_comp);
                Real thv = theta * (1 + 0.61*qv_arr(i,j,k) - ql);
                fab_arr(i, j, k,31) = w_cc_arr(i,j,k) * thv; // w*thv
            });
        } // mfi
    } // use_moisture
 
    h_avg_rho   = sumToLine(mf_out, 0,1,domain,zdir);
    h_avg_th    = sumToLine(mf_out, 1,1,domain,zdir);
    h_avg_ksgs  = sumToLine(mf_out, 2,1,domain,zdir);
    h_avg_Kmv   = sumToLine(mf_out, 3,1,domain,zdir);
    h_avg_Khv   = sumToLine(mf_out, 4,1,domain,zdir);
    h_avg_uu    = sumToLine(mf_out, 5,1,domain,zdir);
    h_avg_uv    = sumToLine(mf_out, 6,1,domain,zdir);
    h_avg_uw    = sumToLine(mf_out, 7,1,domain,zdir);
    h_avg_vv    = sumToLine(mf_out, 8,1,domain,zdir);
    h_avg_vw    = sumToLine(mf_out, 9,1,domain,zdir);
    h_avg_ww    = sumToLine(mf_out,10,1,domain,zdir);
    h_avg_uth   = sumToLine(mf_out,11,1,domain,zdir);
    h_avg_vth   = sumToLine(mf_out,12,1,domain,zdir);
    h_avg_wth   = sumToLine(mf_out,13,1,domain,zdir);
    h_avg_thth  = sumToLine(mf_out,14,1,domain,zdir);
    h_avg_uiuiu = sumToLine(mf_out,15,1,domain,zdir);
    h_avg_uiuiv = sumToLine(mf_out,16,1,domain,zdir);
    h_avg_uiuiw = sumToLine(mf_out,17,1,domain,zdir);
    h_avg_p     = sumToLine(mf_out,18,1,domain,zdir);
    h_avg_pu    = sumToLine(mf_out,19,1,domain,zdir);
    h_avg_pv    = sumToLine(mf_out,20,1,domain,zdir);
    h_avg_pw    = sumToLine(mf_out,21,1,domain,zdir);
    h_avg_qv    = sumToLine(mf_out,22,1,domain,zdir);
    h_avg_qc    = sumToLine(mf_out,23,1,domain,zdir);
    h_avg_qr    = sumToLine(mf_out,24,1,domain,zdir);
    h_avg_wqv   = sumToLine(mf_out,25,1,domain,zdir);
    h_avg_wqc   = sumToLine(mf_out,26,1,domain,zdir);
    h_avg_wqr   = sumToLine(mf_out,27,1,domain,zdir);
    h_avg_qi    = sumToLine(mf_out,28,1,domain,zdir);
    h_avg_qs    = sumToLine(mf_out,29,1,domain,zdir);
    h_avg_qg    = sumToLine(mf_out,30,1,domain,zdir);
    h_avg_wthv  = sumToLine(mf_out,31,1,domain,zdir);
 
    // Divide by the total number of cells we are averaging over
    int h_avg_u_size = static_cast<int>(h_avg_u.size());
    for (int k = 0; k < h_avg_u_size; ++k) {
        h_avg_rho[k] /= area_z;
        h_avg_ksgs[k]  /= area_z;
        h_avg_Kmv[k]   /= area_z;
        h_avg_Khv[k]   /= area_z;
        h_avg_th[k]    /= area_z;
        h_avg_thth[k]  /= area_z;
        h_avg_uu[k]    /= area_z;
        h_avg_uv[k]    /= area_z;
        h_avg_uw[k]    /= area_z;
        h_avg_vv[k]    /= area_z;
        h_avg_vw[k]    /= area_z;
        h_avg_ww[k]    /= area_z;
        h_avg_uth[k]   /= area_z;
        h_avg_vth[k]   /= area_z;
        h_avg_wth[k]   /= area_z;
        h_avg_uiuiu[k] /= area_z;
        h_avg_uiuiv[k] /= area_z;
        h_avg_uiuiw[k] /= area_z;
        h_avg_p[k]     /= area_z;
        h_avg_pu[k]    /= area_z;
        h_avg_pv[k]    /= area_z;
        h_avg_pw[k]    /= area_z;
        h_avg_qv[k]    /= area_z;
        h_avg_qc[k]    /= area_z;
        h_avg_qr[k]    /= area_z;
        h_avg_wqv[k]   /= area_z;
        h_avg_wqc[k]   /= area_z;
        h_avg_wqr[k]   /= area_z;
        h_avg_qi[k]    /= area_z;
        h_avg_qs[k]    /= area_z;
        h_avg_qg[k]    /= area_z;
        h_avg_wthv[k]  /= area_z;
    }
 
#if 0
    // Here we print the integrated total kinetic energy as computed in the 1D profile above
    Real sum = 0.;
    Real dz = geom[0].ProbHi(2) / static_cast<Real>(h_avg_u_size);
    for (int k = 0; k < h_avg_u_size; ++k) {
        sum += h_avg_kturb[k] * h_avg_rho[k] * dz;
    }
    amrex::Print() << "ITKE " << time << " " << sum << " using " << h_avg_u_size << " " << dz << std::endl;
#endif
}

Here is the call graph for this function:

derive_diag_profiles_stag()

void ERF::derive_diag_profiles_stag	(	amrex::Real	time,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_u,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_v,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_w,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_rho,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_th,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_ksgs,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_Kmv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_Khv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qc,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qr,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wqv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wqc,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wqr,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qi,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qs,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_qg,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uu,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_vv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_vw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_ww,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_uth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_vth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_thth,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_ku,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_kv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_kw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_p,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_pu,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_pv,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_pw,
		amrex::Gpu::HostVector< amrex::Real > &	h_avg_wthv
	)

Computes the profiles for diagnostic quantities at staggered heights.

Parameters

h_avg_u	Profile for x-velocity on Host
h_avg_v	Profile for y-velocity on Host
h_avg_w	Profile for z-velocity on Host
h_avg_rho	Profile for density on Host
h_avg_th	Profile for potential temperature on Host
h_avg_ksgs	Profile for Kinetic Energy on Host
h_avg_uu	Profile for x-velocity squared on Host
h_avg_uv	Profile for x-velocity * y-velocity on Host
h_avg_uw	Profile for x-velocity * z-velocity on Host
h_avg_vv	Profile for y-velocity squared on Host
h_avg_vw	Profile for y-velocity * z-velocity on Host
h_avg_ww	Profile for z-velocity squared on Host
h_avg_uth	Profile for x-velocity * potential temperature on Host
h_avg_uiuiu	Profile for u_i*u_i*u triple product on Host
h_avg_uiuiv	Profile for u_i*u_i*v triple product on Host
h_avg_uiuiw	Profile for u_i*u_i*w triple product on Host
h_avg_p	Profile for pressure perturbation on Host
h_avg_pu	Profile for pressure perturbation * x-velocity on Host
h_avg_pv	Profile for pressure perturbation * y-velocity on Host
h_avg_pw	Profile for pressure perturbation * z-velocity on Host

{
    // We assume that this is always called at level 0
    int lev = 0;
 
    bool l_use_kturb = ((solverChoice.turbChoice[lev].les_type != LESType::None) ||
                        (solverChoice.turbChoice[lev].pbl_type != PBLType::None));
    bool l_use_KE   = ( (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) ||
                        (solverChoice.turbChoice[lev].pbl_type == PBLType::MYNN25) );
 
    // Note: "uiui" == u_i*u_i = u*u + v*v + w*w
    // This will hold rho, theta, ksgs, Kmh, Kmv, uu, uv, vv, uth, vth,
    //       indices:   0      1     2    3    4   5   6   7    8    9
    //                thth, uiuiu, uiuiv, p, pu, pv, qv, qc, qr, qi, qs, qg
    //                  10     11     12 13  14  15  16  17  18  19  20  21
    MultiFab mf_out(grids[lev], dmap[lev], 22, 0);
 
    // This will hold uw, vw, ww, wth, uiuiw, pw, wqv, wqc, wqr, wthv
    //       indices:  0   1   2    3      4   5    6    7    8     9
    MultiFab mf_out_stag(convert(grids[lev], IntVect(0,0,1)), dmap[lev], 10, 0);
 
    // This is only used to average u and v; w is not averaged to cell centers
    MultiFab mf_vels(grids[lev], dmap[lev], 2, 0);
 
    MultiFab  u_cc(mf_vels, make_alias, 0, 1); // u at cell centers
    MultiFab  v_cc(mf_vels, make_alias, 1, 1); // v at cell centers
    MultiFab  w_fc(vars_new[lev][Vars::zvel], make_alias, 0, 1); // w at face centers (staggered)
 
    int zdir = 2;
    auto domain = geom[0].Domain();
    Box stag_domain = domain;
    stag_domain.convert(IntVect(0,0,1));
 
    int nvars = vars_new[lev][Vars::cons].nComp();
    MultiFab mf_cons(vars_new[lev][Vars::cons], make_alias, 0, nvars);
 
    MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1);
 
    bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
 
    for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        const Array4<Real>& fab_arr = mf_out.array(mfi);
        const Array4<Real>& fab_arr_stag = mf_out_stag.array(mfi);
        const Array4<Real>& u_arr = vars_new[lev][Vars::xvel].array(mfi);
        const Array4<Real>& v_arr = vars_new[lev][Vars::yvel].array(mfi);
        const Array4<Real>& u_cc_arr =  u_cc.array(mfi);
        const Array4<Real>& v_cc_arr =  v_cc.array(mfi);
        const Array4<Real>& w_fc_arr =  w_fc.array(mfi);
        const Array4<Real>& cons_arr = mf_cons.array(mfi);
        const Array4<Real>&   p0_arr = p_hse.array(mfi);
        const Array4<const Real>& eta_arr = (l_use_kturb) ? eddyDiffs_lev[lev]->const_array(mfi) :
                                                            Array4<const Real>{};
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            u_cc_arr(i,j,k) = 0.5 * (u_arr(i,j,k) + u_arr(i+1,j  ,k));
            v_cc_arr(i,j,k) = 0.5 * (v_arr(i,j,k) + v_arr(i  ,j+1,k));
 
            Real theta = cons_arr(i,j,k,RhoTheta_comp) / cons_arr(i,j,k,Rho_comp);
            fab_arr(i, j, k, 0) = cons_arr(i,j,k,Rho_comp);
            fab_arr(i, j, k, 1) = theta;
            Real ksgs = 0.0;
            if (l_use_KE) {
                ksgs = cons_arr(i,j,k,RhoKE_comp) / cons_arr(i,j,k,Rho_comp);
            }
            fab_arr(i, j, k, 2) = ksgs;
            if (l_use_kturb) {
                fab_arr(i, j, k, 3) = eta_arr(i,j,k,EddyDiff::Mom_v); // Kmv
                fab_arr(i, j, k, 4) = eta_arr(i,j,k,EddyDiff::Theta_v); // Khv
            } else {
                fab_arr(i, j, k, 3) = 0.0;
                fab_arr(i, j, k, 4) = 0.0;
            }
            fab_arr(i, j, k, 5) = u_cc_arr(i,j,k) * u_cc_arr(i,j,k);   // u*u
            fab_arr(i, j, k, 6) = u_cc_arr(i,j,k) * v_cc_arr(i,j,k);   // u*v
            fab_arr(i, j, k, 7) = v_cc_arr(i,j,k) * v_cc_arr(i,j,k);   // v*v
            fab_arr(i, j, k, 8) = u_cc_arr(i,j,k) * theta;             // u*th
            fab_arr(i, j, k, 9) = v_cc_arr(i,j,k) * theta;             // v*th
            fab_arr(i, j, k,10) = theta * theta;                       // th*th
 
            Real wcc = 0.5 * (w_fc_arr(i,j,k) + w_fc_arr(i,j,k+1));
 
            // if the number of fields is changed above, then be sure to update
            // the following def!
            Real uiui = fab_arr(i,j,k,5) + fab_arr(i,j,k,7) + wcc*wcc;
            fab_arr(i, j, k,11) = uiui * u_cc_arr(i,j,k);           // (ui*ui)*u
            fab_arr(i, j, k,12) = uiui * v_cc_arr(i,j,k);           // (ui*ui)*v
 
            if (!use_moisture) {
                Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp));
                p -= p0_arr(i,j,k);
                fab_arr(i, j, k,13) = p;                       // p
                fab_arr(i, j, k,14) = p * u_cc_arr(i,j,k);     // p*u
                fab_arr(i, j, k,15) = p * v_cc_arr(i,j,k);     // p*v
                fab_arr(i, j, k,16) = 0.;  // qv
                fab_arr(i, j, k,17) = 0.;  // qc
                fab_arr(i, j, k,18) = 0.;  // qr
                fab_arr(i, j, k,19) = 0.;  // qi
                fab_arr(i, j, k,20) = 0.;  // qs
                fab_arr(i, j, k,21) = 0.;  // qg
            }
        });
 
        const Box& zbx = mfi.tilebox(IntVect(0,0,1));
        ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            // average to z faces (first to cell centers, then in z)
            Real uface = 0.25 * ( u_arr(i  ,j,k) + u_arr(i  ,j,k-1)
                                + u_arr(i+1,j,k) + u_arr(i+1,j,k-1));
            Real vface = 0.25 * ( v_arr(i,j  ,k) + v_arr(i,j  ,k-1)
                                + v_arr(i,j+1,k) + v_arr(i,j+1,k-1));
            Real theta0 = cons_arr(i,j,k  ,RhoTheta_comp) / cons_arr(i,j,k  ,Rho_comp);
            Real theta1 = cons_arr(i,j,k-1,RhoTheta_comp) / cons_arr(i,j,k-1,Rho_comp);
            Real thface = 0.5*(theta0 + theta1);
            fab_arr_stag(i,j,k,0) =           uface * w_fc_arr(i,j,k); // u*w
            fab_arr_stag(i,j,k,1) =           vface * w_fc_arr(i,j,k); // v*w
            fab_arr_stag(i,j,k,2) = w_fc_arr(i,j,k) * w_fc_arr(i,j,k); // w*w
            fab_arr_stag(i,j,k,3) =          thface * w_fc_arr(i,j,k); // th*w
            Real uiui = uface*uface + vface*vface + fab_arr_stag(i,j,k,2);
            fab_arr_stag(i,j,k,4) = uiui * w_fc_arr(i,j,k); // (ui*ui)*w
            if (!use_moisture) {
                Real p0 = getPgivenRTh(cons_arr(i, j, k  , RhoTheta_comp)) - p0_arr(i,j,k  );
                Real p1 = getPgivenRTh(cons_arr(i, j, k-1, RhoTheta_comp)) - p0_arr(i,j,k-1);
                Real pface = 0.5 * (p0 + p1);
                fab_arr_stag(i,j,k,5) = pface * w_fc_arr(i,j,k);       // p*w
                fab_arr_stag(i,j,k,6) = 0.;  // w*qv
                fab_arr_stag(i,j,k,7) = 0.;  // w*qc
                fab_arr_stag(i,j,k,8) = 0.;  // w*qr
                fab_arr_stag(i,j,k,9) = 0.;  // w*thv
            }
        });
 
    } // mfi
 
    if (use_moisture)
    {
        int n_qstate = micro->Get_Qstate_Size();
 
        for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            const Box& bx = mfi.tilebox();
            const Array4<Real>& fab_arr  = mf_out.array(mfi);
            const Array4<Real>& fab_arr_stag  = mf_out_stag.array(mfi);
            const Array4<Real>& cons_arr = mf_cons.array(mfi);
            const Array4<Real>& u_cc_arr =  u_cc.array(mfi);
            const Array4<Real>& v_cc_arr =  v_cc.array(mfi);
            const Array4<Real>& w_fc_arr =  w_fc.array(mfi);
            const Array4<Real>&   p0_arr = p_hse.array(mfi);
            const Array4<Real>&   qv_arr = qmoist[0][0]->array(mfi); // TODO: Is this written only on lev 0?
 
            int rhoqr_comp = solverChoice.RhoQr_comp;
 
            ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));
 
                p -= p0_arr(i,j,k);
                fab_arr(i, j, k,13) = p;                       // p
                fab_arr(i, j, k,14) = p * u_cc_arr(i,j,k);     // p*u
                fab_arr(i, j, k,15) = p * v_cc_arr(i,j,k);     // p*v
                fab_arr(i, j, k,16) = cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp);  // qv
                fab_arr(i, j, k,17) = cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp);  // qc
                fab_arr(i, j, k,18) = cons_arr(i,j,k,rhoqr_comp) / cons_arr(i,j,k,Rho_comp);  // qr
                if (n_qstate > 3) { // SAM model
                    fab_arr(i, j, k,19) = cons_arr(i,j,k,RhoQ3_comp) / cons_arr(i,j,k,Rho_comp);  // qi
                    fab_arr(i, j, k,20) = cons_arr(i,j,k,RhoQ5_comp) / cons_arr(i,j,k,Rho_comp);  // qs
                    fab_arr(i, j, k,21) = cons_arr(i,j,k,RhoQ6_comp) / cons_arr(i,j,k,Rho_comp);  // qg
                } else {
                    fab_arr(i, j, k,19) = 0.0;  // qi
                    fab_arr(i, j, k,20) = 0.0;  // qs
                    fab_arr(i, j, k,21) = 0.0;  // qg
                }
            });
 
            const Box& zbx = mfi.tilebox(IntVect(0,0,1));
            ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                Real p0 = getPgivenRTh(cons_arr(i, j, k  , RhoTheta_comp), qv_arr(i,j,k  )) - p0_arr(i,j,k  );
                Real p1 = getPgivenRTh(cons_arr(i, j, k-1, RhoTheta_comp), qv_arr(i,j,k-1)) - p0_arr(i,j,k-1);
                Real pface = 0.5 * (p0 + p1);
 
                Real qv0 = cons_arr(i,j,k  ,RhoQ1_comp) / cons_arr(i,j,k  ,Rho_comp);
                Real qv1 = cons_arr(i,j,k-1,RhoQ1_comp) / cons_arr(i,j,k-1,Rho_comp);
                Real qc0 = cons_arr(i,j,k  ,RhoQ2_comp) / cons_arr(i,j,k  ,Rho_comp);
                Real qc1 = cons_arr(i,j,k-1,RhoQ2_comp) / cons_arr(i,j,k-1,Rho_comp);
                Real qr0 = cons_arr(i,j,k  ,RhoQ3_comp) / cons_arr(i,j,k  ,Rho_comp);
                Real qr1 = cons_arr(i,j,k-1,RhoQ3_comp) / cons_arr(i,j,k-1,Rho_comp);
                Real qvface = 0.5 * (qv0 + qv1);
                Real qcface = 0.5 * (qc0 + qc1);
                Real qrface = 0.5 * (qr0 + qr1);
 
                Real theta0 = cons_arr(i,j,k  ,RhoTheta_comp) / cons_arr(i,j,k  ,Rho_comp);
                Real theta1 = cons_arr(i,j,k-1,RhoTheta_comp) / cons_arr(i,j,k-1,Rho_comp);
                Real thface = 0.5*(theta0 + theta1);
                Real ql = qcface + qrface;
                Real thv = thface * (1 + 0.61*qvface - ql);
 
                fab_arr_stag(i,j,k,5) = pface  * w_fc_arr(i,j,k); // p*w
                fab_arr_stag(i,j,k,6) = qvface * w_fc_arr(i,j,k); // w*qv
                fab_arr_stag(i,j,k,7) = qcface * w_fc_arr(i,j,k); // w*qc
                fab_arr_stag(i,j,k,8) = qrface * w_fc_arr(i,j,k); // w*qr
                fab_arr_stag(i,j,k,9) = thv    * w_fc_arr(i,j,k); // w*thv
            });
        } // mfi
    } // use_moisture
 
    // Sum in the horizontal plane
    h_avg_u = sumToLine(u_cc,0,1,     domain,zdir);
    h_avg_v = sumToLine(v_cc,0,1,     domain,zdir);
    h_avg_w = sumToLine(w_fc,0,1,stag_domain,zdir);
 
    h_avg_rho   = sumToLine(mf_out, 0,1,domain,zdir);
    h_avg_th    = sumToLine(mf_out, 1,1,domain,zdir);
    h_avg_ksgs  = sumToLine(mf_out, 2,1,domain,zdir);
    h_avg_Kmv   = sumToLine(mf_out, 3,1,domain,zdir);
    h_avg_Khv   = sumToLine(mf_out, 4,1,domain,zdir);
    h_avg_uu    = sumToLine(mf_out, 5,1,domain,zdir);
    h_avg_uv    = sumToLine(mf_out, 6,1,domain,zdir);
    h_avg_vv    = sumToLine(mf_out, 7,1,domain,zdir);
    h_avg_uth   = sumToLine(mf_out, 8,1,domain,zdir);
    h_avg_vth   = sumToLine(mf_out, 9,1,domain,zdir);
    h_avg_thth  = sumToLine(mf_out,10,1,domain,zdir);
    h_avg_uiuiu = sumToLine(mf_out,11,1,domain,zdir);
    h_avg_uiuiv = sumToLine(mf_out,12,1,domain,zdir);
    h_avg_p     = sumToLine(mf_out,13,1,domain,zdir);
    h_avg_pu    = sumToLine(mf_out,14,1,domain,zdir);
    h_avg_pv    = sumToLine(mf_out,15,1,domain,zdir);
    h_avg_qv    = sumToLine(mf_out,16,1,domain,zdir);
    h_avg_qc    = sumToLine(mf_out,17,1,domain,zdir);
    h_avg_qr    = sumToLine(mf_out,18,1,domain,zdir);
    h_avg_qi    = sumToLine(mf_out,19,1,domain,zdir);
    h_avg_qs    = sumToLine(mf_out,20,1,domain,zdir);
    h_avg_qg    = sumToLine(mf_out,21,1,domain,zdir);
 
    h_avg_uw    = sumToLine(mf_out_stag,0,1,stag_domain,zdir);
    h_avg_vw    = sumToLine(mf_out_stag,1,1,stag_domain,zdir);
    h_avg_ww    = sumToLine(mf_out_stag,2,1,stag_domain,zdir);
    h_avg_wth   = sumToLine(mf_out_stag,3,1,stag_domain,zdir);
    h_avg_uiuiw = sumToLine(mf_out_stag,4,1,stag_domain,zdir);
    h_avg_pw    = sumToLine(mf_out_stag,5,1,stag_domain,zdir);
    h_avg_wqv   = sumToLine(mf_out_stag,6,1,stag_domain,zdir);
    h_avg_wqc   = sumToLine(mf_out_stag,7,1,stag_domain,zdir);
    h_avg_wqr   = sumToLine(mf_out_stag,8,1,stag_domain,zdir);
    h_avg_wthv  = sumToLine(mf_out_stag,9,1,stag_domain,zdir);
 
    // Divide by the total number of cells we are averaging over
    Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
    int unstag_size = h_avg_w.size() - 1; // _un_staggered heights
    for (int k = 0; k < unstag_size; ++k) {
        h_avg_u[k]     /= area_z;
        h_avg_v[k]     /= area_z;
        h_avg_rho[k]   /= area_z;
        h_avg_ksgs[k]  /= area_z;
        h_avg_Kmv[k]   /= area_z;
        h_avg_Khv[k]   /= area_z;
        h_avg_th[k]    /= area_z;
        h_avg_thth[k]  /= area_z;
        h_avg_uu[k]    /= area_z;
        h_avg_uv[k]    /= area_z;
        h_avg_vv[k]    /= area_z;
        h_avg_uth[k]   /= area_z;
        h_avg_vth[k]   /= area_z;
        h_avg_uiuiu[k] /= area_z;
        h_avg_uiuiv[k] /= area_z;
        h_avg_p[k]     /= area_z;
        h_avg_pu[k]    /= area_z;
        h_avg_pv[k]    /= area_z;
        h_avg_qv[k]    /= area_z;
        h_avg_qc[k]    /= area_z;
        h_avg_qr[k]    /= area_z;
        h_avg_qi[k]    /= area_z;
        h_avg_qs[k]    /= area_z;
        h_avg_qg[k]    /= area_z;
    }
 
    for (int k = 0; k < unstag_size+1; ++k) { // staggered heights
        h_avg_w[k]     /= area_z;
        h_avg_uw[k]    /= area_z;
        h_avg_vw[k]    /= area_z;
        h_avg_ww[k]    /= area_z;
        h_avg_wth[k]   /= area_z;
        h_avg_uiuiw[k] /= area_z;
        h_avg_pw[k]    /= area_z;
        h_avg_wqv[k]   /= area_z;
        h_avg_wqc[k]   /= area_z;
        h_avg_wqr[k]   /= area_z;
        h_avg_wthv[k]  /= area_z;
    }
}

Here is the call graph for this function:

derive_stress_profiles()

void ERF::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
	)

{
    int lev = 0;
 
    // This will hold the stress tensor components
    MultiFab mf_out(grids[lev], dmap[lev], 10, 0);
 
    bool l_use_moist   = ( solverChoice.moisture_type != MoistureType::None );
 
    for ( MFIter mfi(mf_out,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        const Array4<Real>& fab_arr = mf_out.array(mfi);
 
        // NOTE: These are from the last RK stage...
        const Array4<const Real>& tau11_arr = Tau11_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau12_arr = Tau12_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau13_arr = Tau13_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau22_arr = Tau22_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau23_arr = Tau23_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau33_arr = Tau33_lev[lev]->const_array(mfi);
 
        // These should be re-calculated during ERF_slow_rhs_post
        // -- just vertical SFS kinematic heat flux for now
        //const Array4<const Real>& hfx1_arr = SFS_hfx1_lev[lev]->const_array(mfi);
        //const Array4<const Real>& hfx2_arr = SFS_hfx2_lev[lev]->const_array(mfi);
        const Array4<const Real>& hfx3_arr = SFS_hfx3_lev[lev]->const_array(mfi);
        const Array4<const Real>& q1fx3_arr = (l_use_moist) ? SFS_q1fx3_lev[lev]->const_array(mfi) :
                                                              Array4<const Real>{};
        const Array4<const Real>& q2fx3_arr = (l_use_moist) ? SFS_q2fx3_lev[lev]->const_array(mfi) :
                                                              Array4<const Real>{};
        const Array4<const Real>& diss_arr = SFS_diss_lev[lev]->const_array(mfi);
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            fab_arr(i, j, k, 0) = tau11_arr(i,j,k);
            fab_arr(i, j, k, 1) = 0.25 * ( tau12_arr(i,j  ,k) + tau12_arr(i+1,j  ,k)
                                         + tau12_arr(i,j+1,k) + tau12_arr(i+1,j+1,k) );
            fab_arr(i, j, k, 2) = 0.25 * ( tau13_arr(i,j,k  ) + tau13_arr(i+1,j,k)
                                         + tau13_arr(i,j,k+1) + tau13_arr(i+1,j,k+1) );
            fab_arr(i, j, k, 3) = tau22_arr(i,j,k);
            fab_arr(i, j, k, 4) = 0.25 * ( tau23_arr(i,j,k  ) + tau23_arr(i,j+1,k)
                                         + tau23_arr(i,j,k+1) + tau23_arr(i,j+1,k+1) );
            fab_arr(i, j, k, 5) = tau33_arr(i,j,k);
            fab_arr(i, j, k, 6) = 0.5 * ( hfx3_arr(i,j,k) + hfx3_arr(i,j,k+1) );
            fab_arr(i, j, k, 7) = (l_use_moist) ? 0.5 * ( q1fx3_arr(i,j,k) + q1fx3_arr(i,j,k+1) ) : 0.0;
            fab_arr(i, j, k, 8) = (l_use_moist) ? 0.5 * ( q2fx3_arr(i,j,k) + q2fx3_arr(i,j,k+1) ) : 0.0;
            fab_arr(i, j, k, 9) =  diss_arr(i,j,k);
        });
    }
 
    int zdir = 2;
    auto domain = geom[0].Domain();
 
    h_avg_tau11 = sumToLine(mf_out,0,1,domain,zdir);
    h_avg_tau12 = sumToLine(mf_out,1,1,domain,zdir);
    h_avg_tau13 = sumToLine(mf_out,2,1,domain,zdir);
    h_avg_tau22 = sumToLine(mf_out,3,1,domain,zdir);
    h_avg_tau23 = sumToLine(mf_out,4,1,domain,zdir);
    h_avg_tau33 = sumToLine(mf_out,5,1,domain,zdir);
    h_avg_hfx3  = sumToLine(mf_out,6,1,domain,zdir);
    h_avg_q1fx3 = sumToLine(mf_out,7,1,domain,zdir);
    h_avg_q2fx3 = sumToLine(mf_out,8,1,domain,zdir);
    h_avg_diss  = sumToLine(mf_out,9,1,domain,zdir);
 
    int ht_size =  h_avg_tau11.size();
 
    // Divide by the total number of cells we are averaging over
    Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
    for (int k = 0; k < ht_size; ++k) {
        h_avg_tau11[k] /= area_z; h_avg_tau12[k] /= area_z; h_avg_tau13[k] /= area_z;
        h_avg_tau22[k] /= area_z; h_avg_tau23[k] /= area_z;
        h_avg_tau33[k] /= area_z;
        h_avg_hfx3[k] /= area_z;
        h_avg_q1fx3[k] /= area_z; h_avg_q2fx3[k] /= area_z;
        h_avg_diss[k] /= area_z;
    }
}

derive_stress_profiles_stag()

void ERF::derive_stress_profiles_stag	(	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
	)

{
    int lev = 0;
 
    // This will hold the stress tensor components
    MultiFab mf_out(grids[lev], dmap[lev], 10, 0);
 
    // This will hold Tau13 and Tau23
    MultiFab mf_out_stag(convert(grids[lev], IntVect(0,0,1)), dmap[lev], 5, 0);
 
    bool l_use_moist   = ( solverChoice.moisture_type != MoistureType::None );
 
    for ( MFIter mfi(mf_out,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        const Array4<Real>& fab_arr = mf_out.array(mfi);
        const Array4<Real>& fab_arr_stag = mf_out_stag.array(mfi);
 
        // NOTE: These are from the last RK stage...
        const Array4<const Real>& tau11_arr = Tau11_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau12_arr = Tau12_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau13_arr = Tau13_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau22_arr = Tau22_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau23_arr = Tau23_lev[lev]->const_array(mfi);
        const Array4<const Real>& tau33_arr = Tau33_lev[lev]->const_array(mfi);
 
        // These should be re-calculated during ERF_slow_rhs_post
        // -- just vertical SFS kinematic heat flux for now
        //const Array4<const Real>& hfx1_arr = SFS_hfx1_lev[lev]->const_array(mfi);
        //const Array4<const Real>& hfx2_arr = SFS_hfx2_lev[lev]->const_array(mfi);
        const Array4<const Real>& hfx3_arr = SFS_hfx3_lev[lev]->const_array(mfi);
        const Array4<const Real>& q1fx3_arr = (l_use_moist) ? SFS_q1fx3_lev[lev]->const_array(mfi) :
                                                              Array4<const Real>{};
        const Array4<const Real>& q2fx3_arr = (l_use_moist) ? SFS_q2fx3_lev[lev]->const_array(mfi) :
                                                              Array4<const Real>{};
        const Array4<const Real>& diss_arr = SFS_diss_lev[lev]->const_array(mfi);
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            fab_arr(i, j, k, 0) = tau11_arr(i,j,k);
            fab_arr(i, j, k, 1) = 0.25 * ( tau12_arr(i,j  ,k) + tau12_arr(i+1,j  ,k)
                                         + tau12_arr(i,j+1,k) + tau12_arr(i+1,j+1,k) );
//          fab_arr(i, j, k, 2) = 0.25 * ( tau13_arr(i,j,k  ) + tau13_arr(i+1,j,k)
//                                       + tau13_arr(i,j,k+1) + tau13_arr(i+1,j,k+1) );
            fab_arr(i, j, k, 3) = tau22_arr(i,j,k);
//          fab_arr(i, j, k, 4) = 0.25 * ( tau23_arr(i,j,k  ) + tau23_arr(i,j+1,k)
//                                       + tau23_arr(i,j,k+1) + tau23_arr(i,j+1,k+1) );
            fab_arr(i, j, k, 5) = tau33_arr(i,j,k);
//          fab_arr(i, j, k, 6) =  hfx3_arr(i,j,k);
//          fab_arr(i, j, k, 7) =  (l_use_moist) ? q1fx3_arr(i,j,k) : 0.0;
//          fab_arr(i, j, k, 8) =  (l_use_moist) ? q2fx3_arr(i,j,k) : 0.0;
            fab_arr(i, j, k, 9) =  diss_arr(i,j,k);
        });
 
        const Box& zbx = mfi.tilebox(IntVect(0,0,1));
        ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            // average from edge to face center
            fab_arr_stag(i,j,k,0) = 0.5*(tau13_arr(i,j,k) + tau13_arr(i+1,j  ,k));
            fab_arr_stag(i,j,k,1) = 0.5*(tau23_arr(i,j,k) + tau23_arr(i  ,j+1,k));
 
            fab_arr_stag(i,j,k,2) =  hfx3_arr(i,j,k);
            fab_arr_stag(i,j,k,3) =  (l_use_moist) ? q1fx3_arr(i,j,k) : 0.0;
            fab_arr_stag(i,j,k,4) =  (l_use_moist) ? q2fx3_arr(i,j,k) : 0.0;
        });
    }
 
    int zdir = 2;
    auto domain = geom[0].Domain();
    Box stag_domain = domain;
    stag_domain.convert(IntVect(0,0,1));
 
    h_avg_tau11 = sumToLine(mf_out,0,1,domain,zdir);
    h_avg_tau12 = sumToLine(mf_out,1,1,domain,zdir);
//  h_avg_tau13 = sumToLine(mf_out,2,1,domain,zdir);
    h_avg_tau22 = sumToLine(mf_out,3,1,domain,zdir);
//  h_avg_tau23 = sumToLine(mf_out,4,1,domain,zdir);
    h_avg_tau33 = sumToLine(mf_out,5,1,domain,zdir);
//  h_avg_hfx3  = sumToLine(mf_out,6,1,domain,zdir);
//  h_avg_q1fx3 = sumToLine(mf_out,7,1,domain,zdir);
//  h_avg_q2fx3 = sumToLine(mf_out,8,1,domain,zdir);
    h_avg_diss  = sumToLine(mf_out,9,1,domain,zdir);
 
    h_avg_tau13 = sumToLine(mf_out_stag,0,1,stag_domain,zdir);
    h_avg_tau23 = sumToLine(mf_out_stag,1,1,stag_domain,zdir);
    h_avg_hfx3  = sumToLine(mf_out_stag,2,1,stag_domain,zdir);
    h_avg_q1fx3 = sumToLine(mf_out_stag,3,1,stag_domain,zdir);
    h_avg_q2fx3 = sumToLine(mf_out_stag,4,1,stag_domain,zdir);
 
    int ht_size =  h_avg_tau11.size(); // _un_staggered
 
    // Divide by the total number of cells we are averaging over
    Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
    for (int k = 0; k < ht_size; ++k) {
        h_avg_tau11[k] /= area_z; h_avg_tau12[k] /= area_z;
        h_avg_tau22[k] /= area_z;
        h_avg_tau33[k] /= area_z;
        h_avg_diss[k] /= area_z;
    }
    for (int k = 0; k < ht_size+1; ++k) { // staggered heights
        h_avg_tau13[k] /= area_z;
        h_avg_tau23[k] /= area_z;
        h_avg_hfx3[k] /= area_z;
        h_avg_q1fx3[k] /= area_z; h_avg_q2fx3[k] /= area_z;
    }
}

derive_upwp()

void ERF::derive_upwp

(

amrex::Vector< amrex::Real > &

h_havg

)

erf_enforce_hse()

void ERF::erf_enforce_hse	(	int	lev,
		amrex::MultiFab &	dens,
		amrex::MultiFab &	pres,
		amrex::MultiFab &	pi,
		amrex::MultiFab &	th,
		std::unique_ptr< amrex::MultiFab > &	z_cc
	)

Enforces hydrostatic equilibrium when using terrain.

Parameters

[in]	lev	Integer specifying the current level
[out]	dens	MultiFab storing base state density
[out]	pres	MultiFab storing base state pressure
[out]	pi	MultiFab storing base state Exner function
[in]	z_cc	Pointer to MultiFab storing cell centered z-coordinates

{
    Real l_gravity = solverChoice.gravity;
    bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None);
 
    const auto geomdata = geom[lev].data();
    const Real dz = geomdata.CellSize(2);
 
    for ( MFIter mfi(dens, TileNoZ()); mfi.isValid(); ++mfi )
    {
        // Create a flat box with same horizontal extent but only one cell in vertical
        const Box& tbz = mfi.nodaltilebox(2);
        int klo = tbz.smallEnd(2);
        int khi = tbz.bigEnd(2);
 
        // Note we only grow by 1 because that is how big z_cc is.
        Box b2d = tbz; // Copy constructor
        b2d.grow(0,1);
        b2d.grow(1,1);
        b2d.setRange(2,0);
 
        // We integrate to the first cell (and below) by using rho in this cell
        // If gravity == 0 this is constant pressure
        // If gravity != 0, hence this is a wall, this gives gp0 = dens[0] * gravity
        // (dens_hse*gravity would also be dens[0]*gravity because we use foextrap for rho at k = -1)
        // Note ng_pres_hse = 1
 
       // We start by assuming pressure on the ground is p_0 (in ERF_Constants.H)
       // Note that gravity is positive
 
        Array4<Real>  rho_arr = dens.array(mfi);
        Array4<Real> pres_arr = pres.array(mfi);
        Array4<Real>   pi_arr =   pi.array(mfi);
        Array4<Real>   th_arr = theta.array(mfi);
        Array4<Real> zcc_arr;
        if (l_use_terrain) {
           zcc_arr = z_cc->array(mfi);
        }
 
        const Real rdOcp = solverChoice.rdOcp;
 
        ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
        {
            // Set value at surface from Newton iteration for rho
            if (klo == 0)
            {
                // Physical height of the terrain at cell center
                Real hz;
                if (l_use_terrain) {
                    hz = zcc_arr(i,j,klo);
                } else {
                    hz = 0.5*dz;
                }
 
                pres_arr(i,j,klo) = p_0 - hz * rho_arr(i,j,klo) * l_gravity;
                  pi_arr(i,j,klo) = getExnergivenP(pres_arr(i,j,klo), rdOcp);
                  th_arr(i,j,klo) =getRhoThetagivenP(pres_arr(i,j,klo)) / rho_arr(i,j,klo);
 
                //
                // Set ghost cell with dz and rho at boundary
                // (We will set the rest of the ghost cells in the boundary condition routine)
                //
                pres_arr(i,j,klo-1) = p_0 + hz * rho_arr(i,j,klo) * l_gravity;
                  pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp);
                  th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1);
 
            } else {
 
                // If level > 0 and klo > 0, we need to use the value of pres_arr(i,j,klo-1) which was
                //    filled from FillPatch-ing it.
                Real dz_loc;
                if (l_use_terrain) {
                    dz_loc = (zcc_arr(i,j,klo) - zcc_arr(i,j,klo-1));
                } else {
                    dz_loc = dz;
                }
 
                Real dens_interp = 0.5*(rho_arr(i,j,klo) + rho_arr(i,j,klo-1));
                pres_arr(i,j,klo) = pres_arr(i,j,klo-1) - dz_loc * dens_interp * l_gravity;
 
                pi_arr(i,j,klo  ) = getExnergivenP(pres_arr(i,j,klo  ), rdOcp);
                th_arr(i,j,klo  ) = getRhoThetagivenP(pres_arr(i,j,klo  )) / rho_arr(i,j,klo  );
 
                pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp);
                th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1);
            }
 
            Real dens_interp;
            if (l_use_terrain) {
                for (int k = klo+1; k <= khi; k++) {
                    Real dz_loc = (zcc_arr(i,j,k) - zcc_arr(i,j,k-1));
                    dens_interp = 0.5*(rho_arr(i,j,k) + rho_arr(i,j,k-1));
                    pres_arr(i,j,k) = pres_arr(i,j,k-1) - dz_loc * dens_interp * l_gravity;
                    pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp);
                    th_arr(i,j,k) = getRhoThetagivenP(pres_arr(i,j,k)) / rho_arr(i,j,k);
                }
            } else {
                for (int k = klo+1; k <= khi; k++) {
                    dens_interp = 0.5*(rho_arr(i,j,k) + rho_arr(i,j,k-1));
                    pres_arr(i,j,k) = pres_arr(i,j,k-1) - dz * dens_interp * l_gravity;
                    pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp);
                    th_arr(i,j,k) = getRhoThetagivenP(pres_arr(i,j,k)) / rho_arr(i,j,k);
                }
            }
        });
 
    } // mfi
 
     dens.FillBoundary(geom[lev].periodicity());
     pres.FillBoundary(geom[lev].periodicity());
       pi.FillBoundary(geom[lev].periodicity());
    theta.FillBoundary(geom[lev].periodicity());
}

Here is the call graph for this function:

ERF_shared()

void ERF::ERF_shared

(

)

{
    if (ParallelDescriptor::IOProcessor()) {
        const char* erf_hash = buildInfoGetGitHash(1);
        const char* amrex_hash = buildInfoGetGitHash(2);
        const char* buildgithash = buildInfoGetBuildGitHash();
        const char* buildgitname = buildInfoGetBuildGitName();
 
        if (strlen(erf_hash) > 0) {
          Print() << "\n"
                         << "ERF git hash: " << erf_hash << "\n";
        }
        if (strlen(amrex_hash) > 0) {
          Print() << "AMReX git hash: " << amrex_hash << "\n";
        }
        if (strlen(buildgithash) > 0) {
          Print() << buildgitname << " git hash: " << buildgithash << "\n";
        }
 
        Print() << "\n";
    }
 
    int nlevs_max = max_level + 1;
 
#ifdef ERF_USE_WINDFARM
    Nturb.resize(nlevs_max);
    vars_windfarm.resize(nlevs_max);
    SMark.resize(nlevs_max);
#endif
 
#if defined(ERF_USE_RRTMGP)
    qheating_rates.resize(nlevs_max);
    sw_lw_fluxes.resize(nlevs_max);
    solar_zenith.resize(nlevs_max);
#endif
 
    // NOTE: size lsm before readparams (chooses the model at all levels)
    lsm.ReSize(nlevs_max);
    lsm_data.resize(nlevs_max);
    lsm_flux.resize(nlevs_max);
 
    // NOTE: size canopy model before readparams (if file exists, we construct)
    m_forest.resize(nlevs_max);
    for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr; }
 
    ReadParameters();
    initializeMicrophysics(nlevs_max);
 
#ifdef ERF_USE_WINDFARM
    initializeWindFarm(nlevs_max);
#endif
 
    const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
    const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
 
    // This is only used when we have flat terrain and stretched z levels.
    stretched_dz_h.resize(nlevs_max);
    stretched_dz_d.resize(nlevs_max);
 
    // Initialize staggered vertical levels for grid stretching or terrain, and
    // to simplify Rayleigh damping layer calculations.
    zlevels_stag.resize(max_level+1);
    init_zlevels(zlevels_stag,
                 stretched_dz_h,
                 stretched_dz_d,
                 geom,
                 refRatio(),
                 solverChoice.grid_stretching_ratio,
                 solverChoice.zsurf,
                 solverChoice.dz0);
 
    if (SolverChoice::terrain_type != TerrainType::None) {
        int nz = geom[0].Domain().length(2) + 1; // staggered
        if (std::fabs(zlevels_stag[0][nz-1]-geom[0].ProbHi(2)) > 1.0e-4) {
            Print() << "Note: prob_hi[2]=" << geom[0].ProbHi(2)
                << " does not match highest requested z level " << zlevels_stag[0][nz-1]
                << std::endl;
        }
        if (std::fabs(zlevels_stag[0][0]-geom[0].ProbLo(2)) > 1.0e-4) {
            Print() << "Note: prob_lo[2]=" << geom[0].ProbLo(2)
                << " does not match lowest requested level " << zlevels_stag[0][0]
                << std::endl;
        }
 
        // Redefine the problem domain here?
    }
 
    prob = amrex_probinit(geom[0].ProbLo(),geom[0].ProbHi());
 
    // Geometry on all levels has been defined already.
 
    // No valid BoxArray and DistributionMapping have been defined.
    // But the arrays for them have been resized.
 
    istep.resize(nlevs_max, 0);
    nsubsteps.resize(nlevs_max, 1);
    for (int lev = 1; lev <= max_level; ++lev) {
        nsubsteps[lev] = MaxRefRatio(lev-1);
    }
 
    t_new.resize(nlevs_max, 0.0);
    t_old.resize(nlevs_max, -1.e100);
    dt.resize(nlevs_max, 1.e100);
    dt_mri_ratio.resize(nlevs_max, 1);
 
    vars_new.resize(nlevs_max);
    vars_old.resize(nlevs_max);
 
    // We resize this regardless in order to pass it without error
    pp_inc.resize(nlevs_max);
 
    rU_new.resize(nlevs_max);
    rV_new.resize(nlevs_max);
    rW_new.resize(nlevs_max);
 
    rU_old.resize(nlevs_max);
    rV_old.resize(nlevs_max);
    rW_old.resize(nlevs_max);
 
    // xmom_crse_rhs.resize(nlevs_max);
    // ymom_crse_rhs.resize(nlevs_max);
    zmom_crse_rhs.resize(nlevs_max);
 
    for (int lev = 0; lev < nlevs_max; ++lev) {
        vars_new[lev].resize(Vars::NumTypes);
        vars_old[lev].resize(Vars::NumTypes);
    }
 
    // Time integrator
    mri_integrator_mem.resize(nlevs_max);
 
    // Physical boundary conditions
    physbcs_cons.resize(nlevs_max);
    physbcs_u.resize(nlevs_max);
    physbcs_v.resize(nlevs_max);
    physbcs_w.resize(nlevs_max);
    physbcs_base.resize(nlevs_max);
 
    // Planes to hold Dirichlet values at boundaries
    xvel_bc_data.resize(nlevs_max);
    yvel_bc_data.resize(nlevs_max);
    zvel_bc_data.resize(nlevs_max);
 
    advflux_reg.resize(nlevs_max);
 
    // Stresses
    Tau11_lev.resize(nlevs_max); Tau22_lev.resize(nlevs_max); Tau33_lev.resize(nlevs_max);
    Tau12_lev.resize(nlevs_max); Tau21_lev.resize(nlevs_max);
    Tau13_lev.resize(nlevs_max); Tau31_lev.resize(nlevs_max);
    Tau23_lev.resize(nlevs_max); Tau32_lev.resize(nlevs_max);
    SFS_hfx1_lev.resize(nlevs_max); SFS_hfx2_lev.resize(nlevs_max); SFS_hfx3_lev.resize(nlevs_max);
    SFS_diss_lev.resize(nlevs_max);
    SFS_q1fx1_lev.resize(nlevs_max); SFS_q1fx2_lev.resize(nlevs_max); SFS_q1fx3_lev.resize(nlevs_max);
    SFS_q2fx3_lev.resize(nlevs_max);
    eddyDiffs_lev.resize(nlevs_max);
    SmnSmn_lev.resize(nlevs_max);
 
    // Sea surface temps
    sst_lev.resize(nlevs_max);
    lmask_lev.resize(nlevs_max);
 
    // Metric terms
    z_phys_nd.resize(nlevs_max);
    z_phys_cc.resize(nlevs_max);
    detJ_cc.resize(nlevs_max);
    ax.resize(nlevs_max);
    ay.resize(nlevs_max);
    az.resize(nlevs_max);
 
    z_phys_nd_new.resize(nlevs_max);
    detJ_cc_new.resize(nlevs_max);
    ax_new.resize(nlevs_max);
    ay_new.resize(nlevs_max);
    az_new.resize(nlevs_max);
 
    z_phys_nd_src.resize(nlevs_max);
    detJ_cc_src.resize(nlevs_max);
    ax_src.resize(nlevs_max);
    ay_src.resize(nlevs_max);
    az_src.resize(nlevs_max);
 
    z_t_rk.resize(nlevs_max);
 
    // Mapfactors
    mapfac_m.resize(nlevs_max);
    mapfac_u.resize(nlevs_max);
    mapfac_v.resize(nlevs_max);
 
    // Thin immersed body
    xflux_imask.resize(nlevs_max);
    yflux_imask.resize(nlevs_max);
    zflux_imask.resize(nlevs_max);
    //overset_imask.resize(nlevs_max);
    thin_xforce.resize(nlevs_max);
    thin_yforce.resize(nlevs_max);
    thin_zforce.resize(nlevs_max);
 
    // Base state
    base_state.resize(nlevs_max);
    base_state_new.resize(nlevs_max);
 
    // Wave coupling data
    Hwave.resize(nlevs_max);
    Lwave.resize(nlevs_max);
    for (int lev = 0; lev < max_level; ++lev)
    {
        Hwave[lev] = nullptr;
        Lwave[lev] = nullptr;
    }
    Hwave_onegrid.resize(nlevs_max);
    Lwave_onegrid.resize(nlevs_max);
    for (int lev = 0; lev < max_level; ++lev)
    {
        Hwave_onegrid[lev] = nullptr;
        Lwave_onegrid[lev] = nullptr;
    }
 
    // Theta prim for MOST
    Theta_prim.resize(nlevs_max);
 
    // Qv prim for MOST
    Qv_prim.resize(nlevs_max);
 
    // Qr prim for MOST
    Qr_prim.resize(nlevs_max);
 
    // Time averaged velocity field
    vel_t_avg.resize(nlevs_max);
    t_avg_cnt.resize(nlevs_max);
 
#ifdef ERF_USE_NETCDF
    // Size lat long arrays if using netcdf
    lat_m.resize(nlevs_max);
    lon_m.resize(nlevs_max);
    for (int lev = 0; lev < max_level; ++lev)
    {
        lat_m[lev] = nullptr;
        lon_m[lev] = nullptr;
    }
#endif
 
    // Initialize tagging criteria for mesh refinement
    refinement_criteria_setup();
 
    for (int lev = 0; lev < max_level; ++lev)
    {
       Print() << "Refinement ratio at level " << lev+1 << " set to be " <<
          ref_ratio[lev][0]  << " " << ref_ratio[lev][1]  <<  " " << ref_ratio[lev][2] << std::endl;
    }
 
    // We will create each of these in MakeNewLevel.../RemakeLevel
    m_factory.resize(max_level+1);
 
#ifdef ERF_USE_EB
    // This is needed before initializing level MultiFabs
    MakeEBGeometry();
#endif
}

Here is the call graph for this function:

ErrorEst()

void ERF::ErrorEst ( int  lev, amrex::TagBoxArray &  tags, amrex::Real  time, int  ngrow  )

override

Function to tag cells for refinement – this overrides the pure virtual function in AmrCore

Parameters

[in]	levc	level of refinement at which we tag cells (0 is coarsest level)
[out]	tags	array of tagged cells
[in]	time	current time

{
    const int clearval = TagBox::CLEAR;
    const int   tagval = TagBox::SET;
 
    for (int j=0; j < ref_tags.size(); ++j)
    {
        std::unique_ptr<MultiFab> mf = std::make_unique<MultiFab>(grids[levc], dmap[levc], 1, 0);
 
        // This allows dynamic refinement based on the value of the density
        if (ref_tags[j].Field() == "density")
        {
            MultiFab::Copy(*mf,vars_new[levc][Vars::cons],Rho_comp,0,1,0);
 
        // This allows dynamic refinement based on the value of qv
        } else if ( ref_tags[j].Field() == "qv" ) {
            MultiFab::Copy(  *mf, vars_new[levc][Vars::cons], RhoQ1_comp, 0, 1, 0);
            MultiFab::Divide(*mf, vars_new[levc][Vars::cons],   Rho_comp, 0, 1, 0);
 
 
        // This allows dynamic refinement based on the value of qc
        } else if (ref_tags[j].Field() == "qc" ) {
            MultiFab::Copy(  *mf, vars_new[levc][Vars::cons], RhoQ2_comp, 0, 1, 0);
            MultiFab::Divide(*mf, vars_new[levc][Vars::cons],   Rho_comp, 0, 1, 0);
 
 
        // This allows dynamic refinement based on the value of the scalar/pressure/theta
        } else if ( (ref_tags[j].Field() == "scalar"  ) ||
                    (ref_tags[j].Field() == "pressure") ||
                    (ref_tags[j].Field() == "theta"   ) )
        {
            for (MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
            {
                const Box& bx = mfi.tilebox();
                auto& dfab = (*mf)[mfi];
                auto& sfab = vars_new[levc][Vars::cons][mfi];
                if (ref_tags[j].Field() == "scalar") {
                    derived::erf_derscalar(bx, dfab, 0, 1, sfab, Geom(levc), time, nullptr, levc);
                } else if (ref_tags[j].Field() == "theta") {
                    derived::erf_dertheta(bx, dfab, 0, 1, sfab, Geom(levc), time, nullptr, levc);
                }
            } // mfi
#ifdef ERF_USE_PARTICLES
        } else {
            //
            // This allows dynamic refinement based on the number of particles per cell
            //
            // Note that we must count all the particles in levels both at and above the current,
            //      since otherwise, e.g., if the particles are all at level 1, counting particles at
            //      level 0 will not trigger refinement when regridding so level 1 will disappear,
            //      then come back at the next regridding
            //
            const auto& particles_namelist( particleData.getNames() );
            mf->setVal(0.0);
            for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++)
            {
                std::string tmp_string(particles_namelist[i]+"_count");
                IntVect rr = IntVect::TheUnitVector();
                if (ref_tags[j].Field() == tmp_string) {
                    for (int lev = levc; lev <= finest_level; lev++)
                    {
                        MultiFab temp_dat(grids[lev], dmap[lev], 1, 0); temp_dat.setVal(0);
                        particleData[particles_namelist[i]]->IncrementWithTotal(temp_dat, lev);
 
                        MultiFab temp_dat_crse(grids[levc], dmap[levc], 1, 0); temp_dat_crse.setVal(0);
 
                        if (lev == levc) {
                            MultiFab::Copy(*mf, temp_dat, 0, 0, 1, 0);
                        } else {
                            for (int d = 0; d < AMREX_SPACEDIM; d++) {
                                rr[d] *= ref_ratio[levc][d];
                            }
                            average_down(temp_dat, temp_dat_crse, 0, 1, rr);
                            MultiFab::Add(*mf, temp_dat_crse, 0, 0, 1, 0);
                        }
                    }
                }
            }
#endif
        }
 
        ref_tags[j](tags,mf.get(),clearval,tagval,time,levc,geom[levc]);
    } // loop over j
}

Here is the call graph for this function:

estTimeStep()

Real ERF::estTimeStep	(	int	level,
		long &	dt_fast_ratio
	)		const

Function that calls estTimeStep for each level

Parameters

[in]	level	level of refinement (coarsest level i 0)
[out]	dt_fast_ratio	ratio of slow to fast time step

{
    BL_PROFILE("ERF::estTimeStep()");
 
    Real estdt_comp = 1.e20;
    Real estdt_lowM = 1.e20;
 
    // We intentionally use the level 0 domain to compute whether to use this direction in the dt calculation
    const int nxc = geom[0].Domain().length(0);
    const int nyc = geom[0].Domain().length(1);
 
    auto const dxinv = geom[level].InvCellSizeArray();
    auto const dzinv = 1.0 / dz_min[level];
 
    MultiFab const& S_new = vars_new[level][Vars::cons];
 
    MultiFab ccvel(grids[level],dmap[level],3,0);
 
    average_face_to_cellcenter(ccvel,0,
                               Array<const MultiFab*,3>{&vars_new[level][Vars::xvel],
                                                        &vars_new[level][Vars::yvel],
                                                        &vars_new[level][Vars::zvel]});
 
    int l_implicit_substepping = (solverChoice.substepping_type[level] == SubsteppingType::Implicit);
    int l_anelastic      = solverChoice.anelastic[level];
 
#ifdef ERF_USE_EB
    EBFArrayBoxFactory ebfact = EBFactory(level);
    const MultiFab& detJ = ebfact.getVolFrac();
#endif
 
#ifdef ERF_USE_EB
    Real estdt_comp_inv = ReduceMax(S_new, ccvel, detJ, 0,
       [=] AMREX_GPU_HOST_DEVICE (Box const& b,
                                  Array4<Real const> const& s,
                                  Array4<Real const> const& u,
                                  Array4<Real const> const& vf) -> Real
#else
    Real estdt_comp_inv = ReduceMax(S_new, ccvel, 0,
       [=] AMREX_GPU_HOST_DEVICE (Box const& b,
                                  Array4<Real const> const& s,
                                  Array4<Real const> const& u) -> Real
#endif
       {
           Real new_comp_dt = -1.e100;
           amrex::Loop(b, [=,&new_comp_dt] (int i, int j, int k) noexcept
           {
#ifdef ERF_USE_EB
               if (vf(i,j,k) > 0.)
#endif
               {
                   const Real rho      = s(i, j, k, Rho_comp);
                   const Real rhotheta = s(i, j, k, RhoTheta_comp);
 
                   // NOTE: even when moisture is present,
                   //       we only use the partial pressure of the dry air
                   //       to compute the soundspeed
                   Real pressure = getPgivenRTh(rhotheta);
                   Real c = std::sqrt(Gamma * pressure / rho);
 
                   // If we are doing implicit acoustic substepping, then the z-direction does not contribute
                   //    to the computation of the time step
                   if (l_implicit_substepping) {
                       if (nxc > 1 && nyc > 1) {
                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
                                                    ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt);
                       } else if (nxc > 1) {
                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), new_comp_dt);
                       } else if (nyc > 1) {
                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt);
                       } else {
                           amrex::Abort("Not sure how to compute dt for this case");
                       }
 
                   // If we are not doing implicit acoustic substepping, then the z-direction contributes
                   //    to the computation of the time step
                   } else {
                       if (nxc > 1 && nyc > 1) {
                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
                                                    ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]),
                                                    ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
                       } else if (nxc > 1) {
                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
                                                    ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
                       } else if (nyc > 1) {
                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]),
                                                    ((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
                       } else {
                           new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,2))+c)*dzinv   ), new_comp_dt);
                       }
 
                   }
               }
           });
           return new_comp_dt;
       });
 
    ParallelDescriptor::ReduceRealMax(estdt_comp_inv);
    estdt_comp = cfl / estdt_comp_inv;
 
     Real estdt_lowM_inv = ReduceMax(ccvel, 0,
       [=] AMREX_GPU_HOST_DEVICE (Box const& b,
                                  Array4<Real const> const& u) -> Real
       {
           Real new_lm_dt = -1.e100;
           Loop(b, [=,&new_lm_dt] (int i, int j, int k) noexcept
           {
               new_lm_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0)))*dxinv[0]),
                                      ((amrex::Math::abs(u(i,j,k,1)))*dxinv[1]),
                                      ((amrex::Math::abs(u(i,j,k,2)))*dxinv[2]), new_lm_dt);
           });
           return new_lm_dt;
       });
 
     ParallelDescriptor::ReduceRealMax(estdt_lowM_inv);
     if (estdt_lowM_inv > 0.0_rt)
         estdt_lowM = cfl / estdt_lowM_inv;
 
     if (verbose) {
         if (fixed_dt[level] <= 0.0) {
             Print() << "Using cfl = " << cfl << " and dx/dy/dz_min = " <<
               1.0/dxinv[0] << " " << 1.0/dxinv[1] << " " << dz_min[level] << std::endl;
             Print() << "Compressible dt at level " << level << ":  " << estdt_comp << std::endl;
             if (estdt_lowM_inv > 0.0_rt) {
                 Print() << "Anelastic dt at level " << level << ":  " << estdt_lowM << std::endl;
             } else {
                 Print() << "Anelastic dt at level " << level << ": undefined " << std::endl;
             }
         }
 
         if (fixed_dt[level] > 0.0) {
             Print() << "Based on cfl of 1.0 " << std::endl;
             Print() << "Compressible dt at level " << level << " would be:  " << estdt_comp/cfl << std::endl;
             if (estdt_lowM_inv > 0.0_rt) {
                 Print() << "Anelastic dt at level " << level << " would be:  " << estdt_lowM/cfl << std::endl;
             } else {
                 Print() << "Anelastic dt at level " << level << " would be undefined " << std::endl;
             }
             Print() << "Fixed dt at level " << level << "       is:  " << fixed_dt[level] << std::endl;
             if (fixed_fast_dt[level] > 0.0) {
                 Print() << "Fixed fast dt at level " << level << "       is:  " << fixed_fast_dt[level] << std::endl;
             }
         }
     }
 
     if (solverChoice.substepping_type[level] != SubsteppingType::None) {
         if (fixed_dt[level] > 0. && fixed_fast_dt[level] > 0.) {
             dt_fast_ratio = static_cast<long>( fixed_dt[level] / fixed_fast_dt[level] );
         } else if (fixed_dt[level] > 0.) {
             dt_fast_ratio = static_cast<long>( 6 );
         }
 
         // Force time step ratio to be an even value
         if (solverChoice.force_stage1_single_substep) {
             if ( dt_fast_ratio%2 != 0) dt_fast_ratio += 1;
         } else {
             if ( dt_fast_ratio%6 != 0) {
                 Print() << "mri_dt_ratio = " << dt_fast_ratio
                         << " not divisible by 6 for N/3 substeps in stage 1" << std::endl;
                 dt_fast_ratio = static_cast<int>(std::ceil(dt_fast_ratio/6.0) * 6);
             }
         }
 
         if (verbose) {
             Print() << "smallest even ratio is: " << dt_fast_ratio << std::endl;
         }
     } // if substepping, either explicit or implicit
 
     if (fixed_dt[level] > 0.0) {
         return fixed_dt[level];
     } else {
         // Anelastic (substepping is not allowed)
         if (l_anelastic) {
             return estdt_lowM;
 
         // Compressible with or without substepping
         } else {
             return estdt_comp;
         }
     }
}

Here is the call graph for this function:

Evolve()

void ERF::Evolve

(

)

{
    BL_PROFILE_VAR("ERF::Evolve()", evolve);
 
    Real cur_time = t_new[0];
 
    // Take one coarse timestep by calling timeStep -- which recursively calls timeStep
    //      for finer levels (with or without subcycling)
    for (int step = istep[0]; step < max_step && cur_time < stop_time; ++step)
    {
        Print() << "\nCoarse STEP " << step+1 << " starts ..." << std::endl;
 
        ComputeDt(step);
 
        // Make sure we have read enough of the boundary plane data to make it through this timestep
        if (input_bndry_planes)
        {
            m_r2d->read_input_files(cur_time,dt[0],m_bc_extdir_vals);
        }
 
        int lev = 0;
        int iteration = 1;
        timeStep(lev, cur_time, iteration);
 
        cur_time  += dt[0];
 
        Print() << "Coarse STEP " << step+1 << " ends." << " TIME = " << cur_time
                << " DT = " << dt[0]  << std::endl;
 
        post_timestep(step, cur_time, dt[0]);
 
        if (writeNow(cur_time, dt[0], step+1, m_plot_int_1, m_plot_per_1)) {
            last_plot_file_step_1 = step+1;
            WritePlotFile(1,plotfile_type_1,plot_var_names_1);
        }
        if (writeNow(cur_time, dt[0], step+1, m_plot_int_2, m_plot_per_2)) {
            last_plot_file_step_2 = step+1;
            WritePlotFile(2,plotfile_type_2,plot_var_names_2);
        }
 
        if (writeNow(cur_time, dt[0], step+1, m_check_int, m_check_per)) {
            last_check_file_step = step+1;
#ifdef ERF_USE_NETCDF
            if (check_type == "netcdf") {
               WriteNCCheckpointFile();
            }
#endif
            if (check_type == "native") {
               WriteCheckpointFile();
            }
        }
 
#ifdef AMREX_MEM_PROFILING
        {
            std::ostringstream ss;
            ss << "[STEP " << step+1 << "]";
            MemProfiler::report(ss.str());
        }
#endif
 
        if (cur_time >= stop_time - 1.e-6*dt[0]) break;
    }
 
    // Write plotfiles at final time
    if ( (m_plot_int_1 > 0 || m_plot_per_1 > 0.) && istep[0] > last_plot_file_step_1 ) {
        WritePlotFile(1,plotfile_type_1,plot_var_names_1);
    }
    if ( (m_plot_int_2 > 0 || m_plot_per_2 > 0.) && istep[0] > last_plot_file_step_2) {
        WritePlotFile(2,plotfile_type_1,plot_var_names_2);
    }
 
    if ( (m_check_int > 0 || m_check_per > 0.) && istep[0] > last_check_file_step) {
#ifdef ERF_USE_NETCDF
        if (check_type == "netcdf") {
           WriteNCCheckpointFile();
        }
#endif
        if (check_type == "native") {
           WriteCheckpointFile();
        }
    }
 
    BL_PROFILE_VAR_STOP(evolve);
}

Referenced by main().

Here is the caller graph for this function:

fill_from_bndryregs()

void ERF::fill_from_bndryregs	(	const amrex::Vector< amrex::MultiFab * > &	mfs,
		amrex::Real	time
	)

{
    //
    // We now assume that if we read in on one face, we read in on all faces
    //
    AMREX_ALWAYS_ASSERT(m_r2d);
 
    int lev = 0;
    const Box& domain = geom[lev].Domain();
 
    const auto& dom_lo = lbound(domain);
    const auto& dom_hi = ubound(domain);
 
    Vector<std::unique_ptr<PlaneVector>>& bndry_data = m_r2d->interp_in_time(time);
 
    const BCRec* bc_ptr = domain_bcs_type_d.data();
 
    // xlo: ori = 0
    // ylo: ori = 1
    // zlo: ori = 2
    // xhi: ori = 3
    // yhi: ori = 4
    // zhi: ori = 5
    const auto& bdatxlo = (*bndry_data[0])[lev].const_array();
    const auto& bdatylo = (*bndry_data[1])[lev].const_array();
    const auto& bdatxhi = (*bndry_data[3])[lev].const_array();
    const auto& bdatyhi = (*bndry_data[4])[lev].const_array();
 
    int bccomp;
 
    for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx)
    {
        MultiFab& mf = *mfs[var_idx];
        const int icomp = 0;
        const int ncomp = mf.nComp();
 
        if (var_idx == Vars::xvel) {
           bccomp = BCVars::xvel_bc;
        } else if (var_idx == Vars::yvel) {
           bccomp = BCVars::yvel_bc;
        } else if (var_idx == Vars::zvel) {
           bccomp = BCVars::zvel_bc;
        } else  if (var_idx == Vars::cons) {
           bccomp = BCVars::cons_bc;
        }
 
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
        for (MFIter mfi(mf); mfi.isValid(); ++mfi)
        {
            const Array4<Real>& dest_arr = mf.array(mfi);
            Box bx = mfi.growntilebox();
 
            // x-faces
            {
            Box bx_xlo(bx); bx_xlo.setBig(0,dom_lo.x-1);
            if (var_idx == Vars::xvel) bx_xlo.setBig(0,dom_lo.x);
 
            Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
            if (var_idx == Vars::xvel) bx_xhi.setSmall(0,dom_hi.x);
 
            ParallelFor(
                bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
                    int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ?
                                   BCVars::RhoScalar_bc_comp : icomp+n;
                    if (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_ingested) {
                        int jb = std::min(std::max(j,dom_lo.y),dom_hi.y);
                        int kb = std::min(std::max(k,dom_lo.z),dom_hi.z);
                        dest_arr(i,j,k,icomp+n) = bdatxlo(dom_lo.x-1,jb,kb,bccomp+n);
                    }
                },
                bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
                    int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ?
                                   BCVars::RhoScalar_bc_comp : icomp+n;
                    if (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_ingested) {
                        int jb = std::min(std::max(j,dom_lo.y),dom_hi.y);
                        int kb = std::min(std::max(k,dom_lo.z),dom_hi.z);
                        dest_arr(i,j,k,icomp+n) = bdatxhi(dom_hi.x+1,jb,kb,bccomp+n);
                    }
                }
            );
            } // x-faces
 
            // y-faces
            {
            Box bx_ylo(bx); bx_ylo.setBig  (1,dom_lo.y-1);
            if (var_idx == Vars::yvel) bx_ylo.setBig(1,dom_lo.y);
 
            Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
            if (var_idx == Vars::yvel) bx_yhi.setSmall(1,dom_hi.y);
 
            ParallelFor(
               bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
                    int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ?
                                   BCVars::RhoScalar_bc_comp : icomp+n;
                    if (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_ingested) {
                        int ib = std::min(std::max(i,dom_lo.x),dom_hi.x);
                        int kb = std::min(std::max(k,dom_lo.z),dom_hi.z);
                        dest_arr(i,j,k,icomp+n) = bdatylo(ib,dom_lo.y-1,kb,bccomp+n);
                    }
                },
                bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
                    int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ?
                                   BCVars::RhoScalar_bc_comp : icomp+n;
                    if (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_ingested) {
                        int ib = std::min(std::max(i,dom_lo.x),dom_hi.x);
                        int kb = std::min(std::max(k,dom_lo.z),dom_hi.z);
                        dest_arr(i,j,k,icomp+n) = bdatyhi(ib,dom_hi.y+1,kb,bccomp+n);
                    }
                }
            );
            } // y-faces
        } // mf
    } // var_idx
}

fill_rhs()

void ERF::fill_rhs ( amrex::MultiFab &  rhs_mf, const amrex::MultiFab &  state_mf, amrex::Real  time, const amrex::Geometry &  geom  )

private

FillBdyCCVels()

void ERF::FillBdyCCVels

(

amrex::Vector< amrex::MultiFab > &

mf_cc_vel

)

{
    // Impose bc's at domain boundaries
    for (int lev = 0; lev <= finest_level; ++lev)
    {
        Box domain(Geom(lev).Domain());
 
        int ihi = domain.bigEnd(0);
        int jhi = domain.bigEnd(1);
        int khi = domain.bigEnd(2);
 
        // Impose periodicity first
        mf_cc_vel[lev].FillBoundary(geom[lev].periodicity());
 
        for (MFIter mfi(mf_cc_vel[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            // Note that we don't fill corners here -- only the cells that share a face
            //      with interior cells -- this is all that is needed to calculate vorticity
            const Box& bx = mfi.tilebox();
            const Array4<Real>& vel_arr = mf_cc_vel[lev].array(mfi);
 
            if (!Geom(lev).isPeriodic(0)) {
                // Low-x side
                if (bx.smallEnd(0) <= domain.smallEnd(0)) {
                    Real mult = (phys_bc_type[0] == ERF_BC::no_slip_wall) ? -1. : 1.;
                    ParallelFor(makeSlab(bx,0,0), [=] AMREX_GPU_DEVICE(int , int j, int k) noexcept
                    {
                        vel_arr(-1,j,k,1) = mult*vel_arr(0,j,k,1); // v
                        vel_arr(-1,j,k,2) = mult*vel_arr(0,j,k,2); // w
                    });
                }
 
                // High-x side
                if (bx.bigEnd(0) >= domain.bigEnd(0)) {
                    Real mult = (phys_bc_type[3] == ERF_BC::no_slip_wall) ? -1. : 1.;
                    ParallelFor(makeSlab(bx,0,0), [=] AMREX_GPU_DEVICE(int , int j, int k) noexcept
                    {
                        vel_arr(ihi+1,j,k,1) = mult*vel_arr(ihi,j,k,1); // v
                        vel_arr(ihi+1,j,k,2) = mult*vel_arr(ihi,j,k,2); // w
                    });
                }
            } // !periodic
 
            if (!Geom(lev).isPeriodic(1)) {
                // Low-y side
                if (bx.smallEnd(1) <= domain.smallEnd(1)) {
                    Real mult = (phys_bc_type[1] == ERF_BC::no_slip_wall) ? -1. : 1.;
                    ParallelFor(makeSlab(bx,1,0), [=] AMREX_GPU_DEVICE(int i, int  , int k) noexcept
                    {
                        vel_arr(i,-1,k,0) = mult*vel_arr(i,0,k,0); // u
                        vel_arr(i,-1,k,2) = mult*vel_arr(i,0,k,2); // w
                    });
                }
 
                // High-y side
                if (bx.bigEnd(1) >= domain.bigEnd(1)) {
                    Real mult = (phys_bc_type[4] == ERF_BC::no_slip_wall) ? -1. : 1.;
                    ParallelFor(makeSlab(bx,1,0), [=] AMREX_GPU_DEVICE(int i, int , int k) noexcept
                    {
                        vel_arr(i,jhi+1,k,0) = mult*vel_arr(i,jhi,k,0); // u
                        vel_arr(i,jhi+1,k,2) = mult*-vel_arr(i,jhi,k,2); // w
                    });
                }
            } // !periodic
 
            if (!Geom(lev).isPeriodic(2)) {
                // Low-z side
                if (bx.smallEnd(2) <= domain.smallEnd(2)) {
                    Real mult = (phys_bc_type[2] == ERF_BC::no_slip_wall) ? -1. : 1.;
                    ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
                    {
                        vel_arr(i,j,-1,0) = mult*vel_arr(i,j,0,0); // u
                        vel_arr(i,j,-1,1) = mult*vel_arr(i,j,0,1); // v
                    });
                }
 
                // High-z side
                if (bx.bigEnd(2) >= domain.bigEnd(2)) {
                    Real mult = (phys_bc_type[5] == ERF_BC::no_slip_wall) ? -1. : 1.;
                    ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
                    {
                        vel_arr(i,j,khi+1,0) = mult*vel_arr(i,j,khi,0); // u
                        vel_arr(i,j,khi+1,1) = mult*vel_arr(i,j,khi,1); // v
                    });
                }
            } // !periodic
        } // MFIter
 
    } // lev
}

FillCoarsePatch()

void ERF::FillCoarsePatch ( int  lev, amrex::Real  time  )

private

{
    BL_PROFILE_VAR("FillCoarsePatch()",FillCoarsePatch);
    AMREX_ASSERT(lev > 0);
 
    //
    //****************************************************************************************************************
    // First fill velocities and density at the COARSE level so we can convert velocity to momenta at the COARSE level
    //****************************************************************************************************************
    //
    bool cons_only = false;
    if (lev == 1) {
        FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel],
                                &vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]},
                                cons_only);
    } else {
        FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel],
                                &vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]},
                               {&vars_new[lev-1][Vars::cons],
                                &rU_new[lev-1], &rV_new[lev-1], &rW_new[lev-1]},
                                base_state[lev-1], base_state[lev-1],
                                false, cons_only);
    }
 
    //
    // ************************************************
    // Convert velocity to momentum at the COARSE level
    // ************************************************
    //
    VelocityToMomentum(vars_new[lev-1][Vars::xvel], IntVect{0},
                       vars_new[lev-1][Vars::yvel], IntVect{0},
                       vars_new[lev-1][Vars::zvel], IntVect{0},
                       vars_new[lev-1][Vars::cons],
                         rU_new[lev-1],
                         rV_new[lev-1],
                         rW_new[lev-1],
                       Geom(lev).Domain(),
                       domain_bcs_type);
    //
    // *****************************************************************
    // Interpolate all cell-centered variables from coarse to fine level
    // *****************************************************************
    //
    Interpolater* mapper_c = &cell_cons_interp;
    Interpolater* mapper_f = &face_cons_linear_interp;
 
    //
    //************************************************************************************************
    // Interpolate cell-centered data from coarse to fine level
    // with InterpFromCoarseLevel which ASSUMES that all ghost cells have already been filled
    // ************************************************************************************************
    IntVect ngvect_cons = vars_new[lev][Vars::cons].nGrowVect();
    int      ncomp_cons = vars_new[lev][Vars::cons].nComp();
 
    InterpFromCoarseLevel(vars_new[lev  ][Vars::cons], ngvect_cons, IntVect(0,0,0),
                          vars_new[lev-1][Vars::cons], 0, 0, ncomp_cons,
                          geom[lev-1], geom[lev],
                          refRatio(lev-1), mapper_c, domain_bcs_type, BCVars::cons_bc);
 
    //
    //************************************************************************************************
    // Interpolate x-momentum from coarse to fine level
    // with InterpFromCoarseLevel which ASSUMES that all ghost cells have already been filled
    // ************************************************************************************************
    //
    InterpFromCoarseLevel(rU_new[lev], IntVect{0}, IntVect{0}, rU_new[lev-1], 0, 0, 1,
                          geom[lev-1], geom[lev],
                          refRatio(lev-1), mapper_f, domain_bcs_type, BCVars::xvel_bc);
 
    //
    //************************************************************************************************
    // Interpolate y-momentum from coarse to fine level
    // with InterpFromCoarseLevel which ASSUMES that all ghost cells have already been filled
    // ************************************************************************************************
    //
    InterpFromCoarseLevel(rV_new[lev], IntVect{0}, IntVect{0}, rV_new[lev-1], 0, 0, 1,
                          geom[lev-1], geom[lev],
                          refRatio(lev-1), mapper_f, domain_bcs_type, BCVars::yvel_bc);
 
    //************************************************************************************************
    // Interpolate z-momentum from coarse to fine level
    // with InterpFromCoarseLevel which ASSUMES that all ghost cells have already been filled
    // ************************************************************************************************
    InterpFromCoarseLevel(rW_new[lev],  IntVect{0}, IntVect{0}, rW_new[lev-1], 0, 0, 1,
                          geom[lev-1], geom[lev],
                          refRatio(lev-1), mapper_f, domain_bcs_type, BCVars::zvel_bc);
    //
    // *********************************************************
    // After interpolation of momentum, convert back to velocity
    // *********************************************************
    //
    for (int which_lev = lev-1; which_lev <= lev; which_lev++)
    {
        MomentumToVelocity(vars_new[which_lev][Vars::xvel],
                           vars_new[which_lev][Vars::yvel],
                           vars_new[which_lev][Vars::zvel],
                           vars_new[which_lev][Vars::cons],
                             rU_new[which_lev],
                             rV_new[which_lev],
                             rW_new[which_lev],
                           Geom(lev).Domain(),
                           domain_bcs_type);
    }
 
    // ***************************************************************************
    // Physical bc's at domain boundary
    // ***************************************************************************
    IntVect ngvect_vels = vars_new[lev][Vars::xvel].nGrowVect();
 
    (*physbcs_cons[lev])(vars_new[lev][Vars::cons],0,ncomp_cons,ngvect_cons,time,BCVars::cons_bc,true);
    (   *physbcs_u[lev])(vars_new[lev][Vars::xvel],0,1         ,ngvect_vels,time,BCVars::xvel_bc,true);
    (   *physbcs_v[lev])(vars_new[lev][Vars::yvel],0,1         ,ngvect_vels,time,BCVars::yvel_bc,true);
    (   *physbcs_w[lev])(vars_new[lev][Vars::zvel],vars_new[lev][Vars::xvel],vars_new[lev][Vars::yvel],
                         ngvect_vels,time,BCVars::zvel_bc,true);
 
    // ***************************************************************************
    // Since lev > 0 here we don't worry about m_r2d or wrfbdy data
    // ***************************************************************************
}

Here is the call graph for this function:

FillIntermediatePatch()

void ERF::FillIntermediatePatch ( int  lev, amrex::Real  time, const amrex::Vector< amrex::MultiFab * > &  mfs_vel, const amrex::Vector< amrex::MultiFab * > &  mfs_mom, int  ng_cons, int  ng_vel, bool  cons_only, int  icomp_cons, int  ncomp_cons, bool  allow_most_bcs = true  )

private

{
    BL_PROFILE_VAR("FillIntermediatePatch()",FillIntermediatePatch);
    Interpolater* mapper;
 
    PhysBCFunctNoOp null_bc;
 
    //
    // ***************************************************************************
    // The first thing we do is interpolate the momenta on the "valid" faces of
    // the fine grids (where the interface is coarse/fine not fine/fine) -- this
    // will not be over-written by interpolation below because the FillPatch
    // operators see these as valid faces.  But we must have these interpolated
    // values in the fine data before we call FillPatchTwoLevels.
    //
    // Also -- note that we might be filling values by interpolation at physical boundaries
    //         here but that's ok because we will overwrite those values when we impose
    //         the physical bc's below
    // ***************************************************************************
    if (lev>0) {
        if (cf_set_width > 0) {
            // We note that mfs_vel[Vars::cons] and mfs_mom[Vars::cons] are in fact the same pointer
            FPr_c[lev-1].FillSet(*mfs_vel[Vars::cons], time, null_bc, domain_bcs_type);
        }
        if ( !cons_only && (cf_set_width >= 0) ) {
            FPr_u[lev-1].FillSet(*mfs_mom[IntVars::xmom], time, null_bc, domain_bcs_type);
            FPr_v[lev-1].FillSet(*mfs_mom[IntVars::ymom], time, null_bc, domain_bcs_type);
            FPr_w[lev-1].FillSet(*mfs_mom[IntVars::zmom], time, null_bc, domain_bcs_type);
        }
    }
 
    // amrex::Print() << "LEVEL " << lev << " CONS ONLY   " << cons_only <<
    //                   " ICOMP NCOMP " << icomp_cons << " " << ncomp_cons << " NGHOST " << ng_cons << std::endl;
 
    if (!cons_only) {
        AMREX_ALWAYS_ASSERT(mfs_mom.size() == IntVars::NumTypes);
        AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes);
    }
 
    // Enforce no penetration for thin immersed body
    if (!cons_only) {
        // Enforce no penetration for thin immersed body
        if (xflux_imask[lev]) {
            ApplyMask(*mfs_mom[IntVars::xmom], *xflux_imask[lev]);
        }
        if (yflux_imask[lev]) {
            ApplyMask(*mfs_mom[IntVars::ymom], *yflux_imask[lev]);
        }
        if (zflux_imask[lev]) {
            ApplyMask(*mfs_mom[IntVars::zmom], *zflux_imask[lev]);
        }
    }
 
    //
    // We now start working on conserved quantities + VELOCITY
    //
    if (lev == 0)
    {
        // We don't do anything here because we will call the physbcs routines below,
        // which calls FillBoundary and fills other domain boundary conditions
        // Physical boundaries will be filled below
 
        if (!cons_only)
        {
            // ***************************************************************************
            // We always come in to this call with updated momenta but we need to create updated velocity
            //    in order to impose the rest of the bc's
            // ***************************************************************************
            // This only fills VALID region of velocity
            MomentumToVelocity(*mfs_vel[Vars::xvel], *mfs_vel[Vars::yvel], *mfs_vel[Vars::zvel],
                               *mfs_vel[Vars::cons],
                               *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom],
                                Geom(lev).Domain(), domain_bcs_type);
        }
    }
    else
    {
        //
        // We must fill a temporary then copy it back so we don't double add/subtract
        //
        MultiFab mf(mfs_vel[Vars::cons]->boxArray(),mfs_vel[Vars::cons]->DistributionMap(),
                    mfs_vel[Vars::cons]->nComp()   ,mfs_vel[Vars::cons]->nGrowVect());
        //
        // Set all components to 1.789e19, then copy just the density from *mfs_vel[Vars::cons]
        //
        mf.setVal(1.789e19);
        MultiFab::Copy(mf,*mfs_vel[Vars::cons],Rho_comp,Rho_comp,1,mf.nGrowVect());
 
        Vector<MultiFab*> fmf = {mfs_vel[Vars::cons],mfs_vel[Vars::cons]};
        Vector<MultiFab*> cmf = {&vars_old[lev-1][Vars::cons], &vars_new[lev-1][Vars::cons]};
        Vector<Real> ctime    = {t_old[lev-1], t_new[lev-1]};
        Vector<Real> ftime    = {time,time};
 
        // Impose physical bc's on fine data (note time and 0 are not used)
        bool do_fb = true; bool do_terrain_adjustment = false;
        (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,
                             do_fb, do_terrain_adjustment);
 
        if ( (icomp_cons+ncomp_cons > 1) && (interpolation_type == StateInterpType::Perturbational) )
        {
            // Divide (rho theta) by rho to get theta
            MultiFab::Divide(*mfs_vel[Vars::cons],*mfs_vel[Vars::cons],Rho_comp,RhoTheta_comp,1,IntVect{0});
 
            // Subtract theta_0 from theta
            MultiFab::Subtract(*mfs_vel[Vars::cons],base_state[lev],BaseState::th0_comp,RhoTheta_comp,1,IntVect{0});
 
            if (!amrex::almostEqual(time,ctime[1])) {
                MultiFab::Divide(vars_old[lev-1][Vars::cons], vars_old[lev-1][Vars::cons],
                                 Rho_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
                MultiFab::Subtract(vars_old[lev-1][Vars::cons], base_state[lev-1],
                                   BaseState::th0_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
            }
            if (!amrex::almostEqual(time,ctime[0])) {
                MultiFab::Divide(vars_new[lev-1][Vars::cons], vars_new[lev-1][Vars::cons],
                                 Rho_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());
                MultiFab::Subtract(vars_new[lev-1][Vars::cons], base_state[lev-1],
                                   BaseState::th0_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());
            }
        }
 
        // Subtract rho_0 from rho before we interpolate -- note we only subtract
        //    on valid region of mf since the ghost cells will be filled below
        if (icomp_cons == 0 && (interpolation_type == StateInterpType::Perturbational))
        {
            MultiFab::Subtract(*mfs_vel[Vars::cons],base_state[lev],BaseState::r0_comp,Rho_comp,1,IntVect{0});
 
            if (!amrex::almostEqual(time,ctime[1])) {
                MultiFab::Subtract(vars_old[lev-1][Vars::cons], base_state[lev-1],
                                   BaseState::r0_comp,Rho_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
            }
            if (!amrex::almostEqual(time,ctime[0])) {
                MultiFab::Subtract(vars_new[lev-1][Vars::cons], base_state[lev-1],
                                   BaseState::r0_comp,Rho_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());
            }
        }
 
        // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled
        mapper = &cell_cons_interp;
        FillPatchTwoLevels(mf, IntVect{ng_cons}, IntVect(0,0,0),
                           time, cmf, ctime, fmf, ftime,
                           icomp_cons, icomp_cons, ncomp_cons, geom[lev-1], geom[lev],
                           refRatio(lev-1), mapper, domain_bcs_type,
                           icomp_cons);
 
        if (icomp_cons == 0 && (interpolation_type == StateInterpType::Perturbational))
        {
            // Restore the coarse values to what they were
            if (!amrex::almostEqual(time,ctime[1])) {
                MultiFab::Add(vars_old[lev-1][Vars::cons], base_state[lev-1],
                              BaseState::r0_comp,Rho_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
            }
            if (!amrex::almostEqual(time,ctime[0])) {
                MultiFab::Add(vars_new[lev-1][Vars::cons], base_state[lev-1],
                              BaseState::r0_comp,Rho_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());
            }
 
            // Set values in the cells outside the domain boundary so that we can do the Add
            //     without worrying about uninitialized values outside the domain -- these
            //     will be filled in the physbcs call
            mf.setDomainBndry(1.234e20,Rho_comp,1,geom[lev]);
 
            // Add rho_0 back to rho after we interpolate -- on all the valid + ghost region
            MultiFab::Add(mf, base_state[lev],BaseState::r0_comp,Rho_comp,1,IntVect{ng_cons});
        }
 
        if ( (icomp_cons+ncomp_cons > 1) && (interpolation_type == StateInterpType::Perturbational) )
        {
            // Add theta_0 to theta
            if (!amrex::almostEqual(time,ctime[1])) {
                MultiFab::Add(vars_old[lev-1][Vars::cons], base_state[lev-1],
                              BaseState::th0_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
                MultiFab::Multiply(vars_old[lev-1][Vars::cons], vars_old[lev-1][Vars::cons],
                                   Rho_comp,RhoTheta_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
            }
            if (!amrex::almostEqual(time,ctime[0])) {
                MultiFab::Add(vars_new[lev-1][Vars::cons], base_state[lev-1],
                              BaseState::th0_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());
                MultiFab::Multiply(vars_new[lev-1][Vars::cons], vars_new[lev-1][Vars::cons],
                                   Rho_comp,RhoTheta_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());
            }
 
            // Multiply theta by rho to get (rho theta)
            MultiFab::Multiply(*mfs_vel[Vars::cons],*mfs_vel[Vars::cons],Rho_comp,RhoTheta_comp,1,IntVect{0});
 
            // Add theta_0 to theta
            MultiFab::Add(*mfs_vel[Vars::cons],base_state[lev],BaseState::th0_comp,RhoTheta_comp,1,IntVect{0});
 
            // Add theta_0 back to theta
            MultiFab::Add(mf,base_state[lev],BaseState::th0_comp,RhoTheta_comp,1,IntVect{ng_cons});
 
            // Multiply (theta) by rho to get (rho theta)
            MultiFab::Multiply(mf,mf,Rho_comp,RhoTheta_comp,1,IntVect{ng_cons});
        }
 
        // Make sure to only copy back the components we worked on
        MultiFab::Copy(*mfs_vel[Vars::cons],mf,icomp_cons,icomp_cons,ncomp_cons,IntVect{ng_cons});
 
        // *****************************************************************************************
 
        if (!cons_only)
        {
            // ***************************************************************************
            // We always come in to this call with updated momenta but we need to create updated velocity
            //    in order to impose the rest of the bc's
            // ***************************************************************************
            // This only fills VALID region of velocity
            MomentumToVelocity(*mfs_vel[Vars::xvel], *mfs_vel[Vars::yvel], *mfs_vel[Vars::zvel],
                               *mfs_vel[Vars::cons],
                               *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom],
                                Geom(lev).Domain(), domain_bcs_type);
 
            mapper = &face_cons_linear_interp;
 
            //
            // NOTE: All interpolation here happens on velocities not momenta;
            //       note we only do the interpolation and FillBoundary here,
            //       physical bc's are imposed later
            //
            // NOTE: This will only fill velocity from coarse grid *outside* the fine grids
            //       unlike the FillSet calls above which filled momenta on the coarse/fine bdy
            //
 
            MultiFab& mfu = *mfs_vel[Vars::xvel];
 
            fmf = {&mfu,&mfu};
            cmf = {&vars_old[lev-1][Vars::xvel], &vars_new[lev-1][Vars::xvel]};
 
            // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled
            FillPatchTwoLevels(mfu, IntVect{ng_vel}, IntVect(0,0,0),
                               time, cmf, ctime, fmf, ftime,
                               0, 0, 1, geom[lev-1], geom[lev],
                               refRatio(lev-1), mapper, domain_bcs_type,
                               BCVars::xvel_bc);
 
            // *****************************************************************************************
 
            MultiFab& mfv = *mfs_vel[Vars::yvel];
 
            fmf = {&mfv,&mfv};
            cmf = {&vars_old[lev-1][Vars::yvel], &vars_new[lev-1][Vars::yvel]};
 
            // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled
            FillPatchTwoLevels(mfv, IntVect{ng_vel}, IntVect(0,0,0),
                               time, cmf, ctime, fmf, ftime,
                               0, 0, 1, geom[lev-1], geom[lev],
                               refRatio(lev-1), mapper, domain_bcs_type,
                               BCVars::yvel_bc);
 
            // *****************************************************************************************
 
            MultiFab& mfw = *mfs_vel[Vars::zvel];
 
            fmf = {&mfw,&mfw};
            cmf = {&vars_old[lev-1][Vars::zvel], &vars_new[lev-1][Vars::zvel]};
 
            // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled
            FillPatchTwoLevels(mfw, IntVect{ng_vel}, IntVect(0,0,0),
                               time, cmf, ctime, fmf, ftime,
                               0, 0, 1, geom[lev-1], geom[lev],
                               refRatio(lev-1), mapper, domain_bcs_type,
                               BCVars::zvel_bc);
        } // !cons_only
    } // lev > 0
 
    // ***************************************************************************
    // Physical bc's at domain boundary
    // ***************************************************************************
    IntVect ngvect_cons = IntVect(ng_cons,ng_cons,ng_cons);
    IntVect ngvect_vels = IntVect(ng_vel ,ng_vel ,ng_vel);
 
    bool do_fb = true;
 
#ifdef ERF_USE_NETCDF
    // We call this here because it is an ERF routine
    if (use_real_bcs && (lev==0)) {
        fill_from_realbdy(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels);
        do_fb = false;
    }
#endif
 
    if (m_r2d) fill_from_bndryregs(mfs_vel,time);
 
    // We call this even if init_type == InitType::Real because this routine will fill the vertical bcs
    (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,ngvect_cons,time,BCVars::cons_bc, do_fb);
    if (!cons_only) {
        (*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc, do_fb);
        (*physbcs_v[lev])(*mfs_vel[Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc, do_fb);
        (*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel],
                          ngvect_vels,time,BCVars::zvel_bc, do_fb);
    }
    // ***************************************************************************
 
    // MOST boundary conditions
    if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs) {
        m_most->impose_most_bcs(lev,mfs_vel,
                                Tau11_lev[lev].get(),
                                Tau22_lev[lev].get(),
                                Tau33_lev[lev].get(),
                                Tau12_lev[lev].get(), Tau21_lev[lev].get(),
                                Tau13_lev[lev].get(), Tau31_lev[lev].get(),
                                Tau23_lev[lev].get(), Tau32_lev[lev].get(),
                                SFS_hfx1_lev[lev].get(),
                                SFS_hfx2_lev[lev].get(),
                                SFS_hfx3_lev[lev].get(),
                                SFS_q1fx1_lev[lev].get(),
                                SFS_q1fx2_lev[lev].get(),
                                SFS_q1fx3_lev[lev].get(),
                                z_phys_nd[lev].get());
    }
 
    // We always come in to this call with momenta so we need to leave with momenta!
    // We need to make sure we convert back on all ghost cells/faces because this is
    // how velocity from fine-fine copies (as well as physical and interpolated bcs) will be filled
    if (!cons_only)
    {
        IntVect ngu = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : mfs_vel[Vars::xvel]->nGrowVect();
        IntVect ngv = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : mfs_vel[Vars::yvel]->nGrowVect();
        IntVect ngw = (!solverChoice.use_NumDiff) ? IntVect(1,1,0) : mfs_vel[Vars::zvel]->nGrowVect();
 
        VelocityToMomentum(*mfs_vel[Vars::xvel], ngu,
                           *mfs_vel[Vars::yvel], ngv,
                           *mfs_vel[Vars::zvel], ngw,
                           *mfs_vel[Vars::cons],
                           *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom],
                           Geom(lev).Domain(),
                           domain_bcs_type);
    }
 
    mfs_mom[IntVars::cons]->FillBoundary(geom[lev].periodicity());
    mfs_mom[IntVars::xmom]->FillBoundary(geom[lev].periodicity());
    mfs_mom[IntVars::ymom]->FillBoundary(geom[lev].periodicity());
    mfs_mom[IntVars::zmom]->FillBoundary(geom[lev].periodicity());
}

Here is the call graph for this function:

FillPatch() [1/2]

void ERF::FillPatch ( int  lev, amrex::Real  time, const amrex::Vector< amrex::MultiFab * > &  mfs_vel, bool  cons_only = false  )

private

FillPatch() [2/2]

void ERF::FillPatch ( int  lev, amrex::Real  time, const amrex::Vector< amrex::MultiFab * > &  mfs_vel, const amrex::Vector< amrex::MultiFab * > &  mfs_mom, const amrex::MultiFab &  old_base_state, const amrex::MultiFab &  new_base_state, bool  fillset = true, bool  cons_only = false  )

private

get_projection_bc()

Array< LinOpBCType, AMREX_SPACEDIM > ERF::get_projection_bc ( amrex::Orientation::Side  side ) const

noexcept

{
    amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> r;
    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
        if (geom[0].isPeriodic(dir)) {
            r[dir] = LinOpBCType::Periodic;
        } else {
            auto bc_type = domain_bc_type[Orientation(dir,side)];
            if (bc_type == "Outflow") {
                r[dir] = LinOpBCType::Dirichlet;
            } else
            {
                r[dir] = LinOpBCType::Neumann;
            }
        }
    }
    return r;
}

getAdvFluxReg()

AMREX_FORCE_INLINE amrex::YAFluxRegister* ERF::getAdvFluxReg ( int  lev )

inlineprivate

    {
        return advflux_reg[lev];
    }

getCPUTime()

static amrex::Real ERF::getCPUTime ( )

inlinestaticprivate

    {
      int numCores = amrex::ParallelDescriptor::NProcs();
#ifdef _OPENMP
      numCores = numCores * omp_get_max_threads();
#endif
 
      amrex::Real T =
        numCores * (amrex::ParallelDescriptor::second() - startCPUTime) +
        previousCPUTimeUsed;
 
      return T;
    }

GotoNextLine()

void ERF::GotoNextLine ( std::istream &  is )

staticprivate

Utility to skip to next line in Header file input stream.

{
    constexpr std::streamsize bl_ignore_max { 100000 };
    is.ignore(bl_ignore_max, '\n');
}

init1DArrays()

void ERF::init1DArrays ( )

private

init_bcs()

void ERF::init_bcs ( )

private

Initializes data structures in the ERF class that specify which boundary conditions we are implementing on each face of the domain.

This function also maps the selected boundary condition types (e.g. Outflow, Inflow, Periodic, Dirichlet, ...) to the specific implementation needed for each variable.

Stores this information in both host and device vectors so it is available for GPU kernels.

{
    bool rho_read = false;
    Vector<Real> cons_dir_init(NBCVAR_max,0.0);
    cons_dir_init[BCVars::Rho_bc_comp] = 1.0;
    cons_dir_init[BCVars::RhoTheta_bc_comp] = -1.0;
    auto f = [this,&rho_read] (std::string const& bcid, Orientation ori)
    {
        // These are simply defaults for Dirichlet faces -- they should be over-written below
        m_bc_extdir_vals[BCVars::Rho_bc_comp][ori]       =  1.0;
        m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori]  = -1.0; // It is important to set this negative
                                                                 // because the sign is tested on below
        m_bc_extdir_vals[BCVars::RhoKE_bc_comp][ori]     = 0.0;
        m_bc_extdir_vals[BCVars::RhoScalar_bc_comp][ori] = 0.0;
        m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori]     = 0.0;
        m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori]     = 0.0;
        m_bc_extdir_vals[BCVars::RhoQ3_bc_comp][ori]     = 0.0;
        m_bc_extdir_vals[BCVars::RhoQ4_bc_comp][ori]     = 0.0;
        m_bc_extdir_vals[BCVars::RhoQ5_bc_comp][ori]     = 0.0;
        m_bc_extdir_vals[BCVars::RhoQ6_bc_comp][ori]     = 0.0;
 
        m_bc_extdir_vals[BCVars::xvel_bc][ori] = 0.0; // default
        m_bc_extdir_vals[BCVars::yvel_bc][ori] = 0.0;
        m_bc_extdir_vals[BCVars::zvel_bc][ori] = 0.0;
 
        // These are simply defaults for Neumann gradients -- they should be over-written below
        m_bc_neumann_vals[BCVars::Rho_bc_comp][ori]       = 0.0;
        m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]  = 0.0;
 
        m_bc_neumann_vals[BCVars::RhoKE_bc_comp][ori]     = 0.0;
        m_bc_neumann_vals[BCVars::RhoScalar_bc_comp][ori] = 0.0;
        m_bc_neumann_vals[BCVars::RhoQ1_bc_comp][ori]     = 0.0;
        m_bc_neumann_vals[BCVars::RhoQ2_bc_comp][ori]     = 0.0;
        m_bc_neumann_vals[BCVars::RhoQ3_bc_comp][ori]     = 0.0;
        m_bc_neumann_vals[BCVars::RhoQ4_bc_comp][ori]     = 0.0;
        m_bc_neumann_vals[BCVars::RhoQ5_bc_comp][ori]     = 0.0;
        m_bc_neumann_vals[BCVars::RhoQ6_bc_comp][ori]     = 0.0;
 
        m_bc_neumann_vals[BCVars::xvel_bc][ori] = 0.0;
        m_bc_neumann_vals[BCVars::yvel_bc][ori] = 0.0;
        m_bc_neumann_vals[BCVars::zvel_bc][ori] = 0.0;
 
        std::string pp_text;
        if (pp_prefix == "erf") {
          pp_text = bcid;
        } else {
          pp_text = pp_prefix + "." + bcid;
        }
        ParmParse pp(pp_text);
        std::string bc_type_in = "null";
        pp.query("type", bc_type_in);
        std::string bc_type = amrex::toLower(bc_type_in);
 
        if (bc_type == "symmetry")
        {
            // Print() << bcid << " set to symmetry.\n";
              phys_bc_type[ori] = ERF_BC::symmetry;
            domain_bc_type[ori] = "Symmetry";
        }
        else if (bc_type == "outflow")
        {
            // Print() << bcid << " set to outflow.\n";
              phys_bc_type[ori] = ERF_BC::outflow;
            domain_bc_type[ori] = "Outflow";
        }
        else if (bc_type == "open")
        {
            // Print() << bcid << " set to open.\n";
            AMREX_ASSERT_WITH_MESSAGE((ori.coordDir() != 2), "Open boundary not valid on zlo or zhi!");
              phys_bc_type[ori] = ERF_BC::open;
            domain_bc_type[ori] = "Open";
        }
        else if (bc_type == "ho_outflow")
        {
            phys_bc_type[ori] = ERF_BC::ho_outflow;
            domain_bc_type[ori] = "HO_Outflow";
        }
 
        else if (bc_type == "inflow")
        {
            // Print() << bcid << " set to inflow.\n";
              phys_bc_type[ori] = ERF_BC::inflow;
            domain_bc_type[ori] = "Inflow";
 
            std::vector<Real> v;
            if (input_bndry_planes && m_r2d->ingested_velocity()) {
                m_bc_extdir_vals[BCVars::xvel_bc][ori] = 0.;
                m_bc_extdir_vals[BCVars::yvel_bc][ori] = 0.;
                m_bc_extdir_vals[BCVars::zvel_bc][ori] = 0.;
            } else {
                // Test for input data file if at xlo face
                std::string dirichlet_file;
                auto file_exists = pp.query("dirichlet_file", dirichlet_file);
                if (file_exists) {
                    init_Dirichlet_bc_data(dirichlet_file);
                } else {
                    pp.getarr("velocity", v, 0, AMREX_SPACEDIM);
                    m_bc_extdir_vals[BCVars::xvel_bc][ori] = v[0];
                    m_bc_extdir_vals[BCVars::yvel_bc][ori] = v[1];
                    m_bc_extdir_vals[BCVars::zvel_bc][ori] = v[2];
                }
            }
 
            Real rho_in = 0.;
            if (input_bndry_planes && m_r2d->ingested_density()) {
                m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] = 0.;
            } else {
                if (!pp.query("density", rho_in)) {
                    amrex::Print() << "Using interior values to set conserved vars" << std::endl;
                }
                m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] = rho_in;
            }
            Real theta_in = 0.;
            if (input_bndry_planes && m_r2d->ingested_theta()) {
                m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] = 0.;
            } else {
                if (rho_in > 0) {
                    pp.get("theta", theta_in);
                }
                m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] = rho_in*theta_in;
            }
            Real scalar_in = 0.;
            if (input_bndry_planes && m_r2d->ingested_scalar()) {
                m_bc_extdir_vals[BCVars::RhoScalar_bc_comp][ori] = 0.;
            } else {
                if (pp.query("scalar", scalar_in))
                m_bc_extdir_vals[BCVars::RhoScalar_bc_comp][ori] = rho_in*scalar_in;
            }
 
            if (solverChoice.moisture_type != MoistureType::None) {
                Real qv_in = 0.;
                if (input_bndry_planes && m_r2d->ingested_q1()) {
                    m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = 0.;
                } else {
                    if (pp.query("qv", qv_in))
                    m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = rho_in*qv_in;
                }
                Real qc_in = 0.;
                if (input_bndry_planes && m_r2d->ingested_q2()) {
                    m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = 0.;
                } else {
                    if (pp.query("qc", qc_in))
                    m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = rho_in*qc_in;
                }
            }
 
            Real KE_in = 0.;
            if (input_bndry_planes && m_r2d->ingested_KE()) {
                m_bc_extdir_vals[BCVars::RhoKE_bc_comp][ori] = 0.;
            } else {
                if (pp.query("KE", KE_in))
                m_bc_extdir_vals[BCVars::RhoKE_bc_comp][ori] = rho_in*KE_in;
            }
        }
        else if (bc_type == "noslipwall")
        {
            // Print() << bcid <<" set to no-slip wall.\n";
              phys_bc_type[ori] = ERF_BC::no_slip_wall;
            domain_bc_type[ori] = "NoSlipWall";
 
            std::vector<Real> v;
 
            // The values of m_bc_extdir_vals default to 0.
            // But if we find "velocity" in the inputs file, use those values instead.
            if (pp.queryarr("velocity", v, 0, AMREX_SPACEDIM))
            {
                v[ori.coordDir()] = 0.0;
                m_bc_extdir_vals[BCVars::xvel_bc][ori] = v[0];
                m_bc_extdir_vals[BCVars::yvel_bc][ori] = v[1];
                m_bc_extdir_vals[BCVars::zvel_bc][ori] = v[2];
            }
 
            Real rho_in;
            rho_read = pp.query("density", rho_in);
            if (rho_read)
            {
                m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] = rho_in;
            }
 
            Real theta_in;
            if (pp.query("theta", theta_in))
            {
               m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] = theta_in*m_bc_extdir_vals[BCVars::Rho_bc_comp][ori];
            }
 
            Real theta_grad_in;
            if (pp.query("theta_grad", theta_grad_in))
            {
                m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori] = theta_grad_in;
            }
 
            Real qv_in;
            if (pp.query("qv", qv_in))
            {
                m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = qv_in*m_bc_extdir_vals[BCVars::Rho_bc_comp][ori];
            }
        }
        else if (bc_type == "slipwall")
        {
            // Print() << bcid <<" set to slip wall.\n";
 
              phys_bc_type[ori] = ERF_BC::slip_wall;
            domain_bc_type[ori] = "SlipWall";
 
            Real rho_in;
            if (pp.query("density", rho_in))
            {
                m_bc_extdir_vals[BCVars::Rho_bc_comp][ori] = rho_in;
            }
 
            Real theta_in;
            if (pp.query("theta", theta_in))
            {
               m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] = theta_in*m_bc_extdir_vals[BCVars::Rho_bc_comp][ori];
            }
 
            Real rho_grad_in;
            if (pp.query("density_grad", rho_grad_in))
            {
               m_bc_neumann_vals[BCVars::Rho_bc_comp][ori] = rho_grad_in;
            }
 
            Real theta_grad_in;
            if (pp.query("theta_grad", theta_grad_in))
            {
               m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori] = theta_grad_in;
            }
        }
        else if (bc_type == "most")
        {
              phys_bc_type[ori] = ERF_BC::MOST;
            domain_bc_type[ori] = "MOST";
        }
        else
        {
            phys_bc_type[ori] = ERF_BC::undefined;
        }
 
        if (geom[0].isPeriodic(ori.coordDir())) {
            domain_bc_type[ori] = "Periodic";
            if (phys_bc_type[ori] == ERF_BC::undefined)
            {
                phys_bc_type[ori] = ERF_BC::periodic;
            } else {
                Abort("Wrong BC type for periodic boundary");
            }
        }
 
        if (phys_bc_type[ori] == ERF_BC::undefined)
        {
             Print() << "BC Type specified for face " << bcid << " is " << bc_type_in << std::endl;
             Abort("This BC type is unknown");
        }
    };
 
    f("xlo", Orientation(Direction::x,Orientation::low));
    f("xhi", Orientation(Direction::x,Orientation::high));
    f("ylo", Orientation(Direction::y,Orientation::low));
    f("yhi", Orientation(Direction::y,Orientation::high));
    f("zlo", Orientation(Direction::z,Orientation::low));
    f("zhi", Orientation(Direction::z,Orientation::high));
 
    // *****************************************************************************
    //
    // Here we translate the physical boundary conditions -- one type per face --
    //     into logical boundary conditions for each velocity component
    //
    // *****************************************************************************
    {
        domain_bcs_type.resize(AMREX_SPACEDIM+NBCVAR_max);
        domain_bcs_type_d.resize(AMREX_SPACEDIM+NBCVAR_max);
 
        for (OrientationIter oit; oit; ++oit) {
            Orientation ori = oit();
            int dir = ori.coordDir();
            Orientation::Side side = ori.faceDir();
            auto const bct = phys_bc_type[ori];
            if ( bct == ERF_BC::symmetry )
            {
                if (side == Orientation::low) {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::reflect_even);
                    }
                    domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::reflect_odd);
                } else {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::reflect_even);
                    }
                    domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::reflect_odd);
                }
            }
            else if (bct == ERF_BC::outflow or bct == ERF_BC::ho_outflow )
            {
                if (side == Orientation::low) {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::foextrap);
                    }
                    domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::neumann_int);
                } else {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::foextrap);
                    }
                    domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::neumann_int);
                }
            }
            else if (bct == ERF_BC::open)
            {
                if (side == Orientation::low) {
                    for (int i = 0; i < AMREX_SPACEDIM; i++)
                        domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::open);
                } else {
                    for (int i = 0; i < AMREX_SPACEDIM; i++)
                        domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::open);
                }
            }
            else if (bct == ERF_BC::inflow)
            {
                if (side == Orientation::low) {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir);
                        if (input_bndry_planes && dir < 2 && m_r2d->ingested_velocity()) {
                            domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir_ingested);
                        }
                    }
                } else {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir);
                        if (input_bndry_planes && dir < 2 && m_r2d->ingested_velocity()) {
                            domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir_ingested);
                        }
                    }
                }
            }
            else if (bct == ERF_BC::no_slip_wall)
            {
                if (side == Orientation::low) {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::ext_dir);
                    }
                } else {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::ext_dir);
                    }
                }
            }
            else if (bct == ERF_BC::slip_wall)
            {
                if (side == Orientation::low) {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::foextrap);
                    }
                    // Only normal direction has ext_dir
                    domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::ext_dir);
 
                } else {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::foextrap);
                    }
                    // Only normal direction has ext_dir
                    domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::ext_dir);
                }
            }
            else if (bct == ERF_BC::periodic)
            {
                if (side == Orientation::low) {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::int_dir);
                    }
                } else {
                    for (int i = 0; i < AMREX_SPACEDIM; i++) {
                        domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::int_dir);
                    }
                }
            }
            else if ( bct ==