diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 4b0dc4cee..102e22fad 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -174,7 +174,8 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Utils/MomentumToVelocity.cpp ${SRC_DIR}/Utils/TerrainMetrics.cpp ${SRC_DIR}/Utils/VelocityToMomentum.cpp - ${SRC_DIR}/Utils/InteriorGhostCells.cpp + ${SRC_DIR}/Utils/InteriorGhostCells.cpp + ${SRC_DIR}/Utils/Time_Avg_Vel.cpp ${SRC_DIR}/Microphysics/SAM/Init_SAM.cpp ${SRC_DIR}/Microphysics/SAM/Cloud_SAM.cpp ${SRC_DIR}/Microphysics/SAM/IceFall.cpp diff --git a/Docs/sphinx_doc/Plotfiles.rst b/Docs/sphinx_doc/Plotfiles.rst index 4ff201b9c..146249f4d 100644 --- a/Docs/sphinx_doc/Plotfiles.rst +++ b/Docs/sphinx_doc/Plotfiles.rst @@ -108,7 +108,9 @@ PlotFile Outputs ================ Plotfiles can include the quantities of several simulation parameters as output. -They are summarized in the list below. +They are summarized in the list below. Note that temporally averaged quantities +(e.g., ``u_t_avg, v_t_avg, w_t_avg, umag_t_avg``) require the user to enable the +storage of the time averaged variables with ``erf.time_avg_vel = true``. Output Options -------------- @@ -201,18 +203,34 @@ Output Options | | | | | | +-----------------------------+------------------+ -| **vorticity_x* | x-component of | +| **vorticity_x** | x-component of | | | vorticity | | | | +-----------------------------+------------------+ -| **vorticity_y* | y-component of | +| **vorticity_y** | y-component of | | | vorticity | | | | +-----------------------------+------------------+ -| **vorticity_z* | z-component of | +| **vorticity_z** | z-component of | | | vorticity | | | | +-----------------------------+------------------+ +| **u_t_avg** | time average of | +| | x-component of | +| | velocity | ++-----------------------------+------------------+ +| **v_t_avg** | time average of | +| | y-component of | +| | velocity | ++-----------------------------+------------------+ +| **w_t_avg** | time average of | +| | z-component of | +| | velocity | ++-----------------------------+------------------+ +| **umag_t_avg** | time average of | +| | velocity mag | +| | | ++-----------------------------+------------------+ | **rhoadv_0** | Conserved scalar | | | | | | | diff --git a/Exec/ABL/inputs_vel_avg b/Exec/ABL/inputs_vel_avg new file mode 100644 index 000000000..9b114b86d --- /dev/null +++ b/Exec/ABL/inputs_vel_avg @@ -0,0 +1,65 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 4000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 1024 1024 1024 +amr.n_cell = 64 64 64 + +geometry.is_periodic = 1 1 0 + +zlo.type = "NoSlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.1 # fixed time step depending on grid resolution + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.time_avg_vel = true +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 10 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoKE rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta u_t_avg v_t_avg w_t_avg umag_t_avg + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "None" +erf.les_type = "Deardorff" +erf.Ck = 0.1 +erf.sigma_k = 1.0 +erf.Ce = 0.1 + +erf.init_type = "uniform" +erf.KE_0 = 0.1 # for Deardorff + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + +# Higher values of perturbations lead to instability +# Instability seems to be coming from BC +prob.U_0_Pert_Mag = 0.08 +prob.V_0_Pert_Mag = 0.08 # +prob.W_0_Pert_Mag = 0.0 diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index aba5ecb26..df38dbb62 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -278,6 +278,9 @@ struct SolverChoice { else { amrex::Abort("Dont know this windfarm_type. Currently the only option is Fitch."); } + + // Test if time averaged data is to be output + pp.query("time_avg_vel",time_avg_vel); } void display() @@ -423,6 +426,9 @@ struct SolverChoice { bool custom_geostrophic_profile = false; bool custom_forcing_prim_vars = false; + // User wishes to output time averaged velocity fields + bool time_avg_vel = false; + // Numerical diffusion bool use_NumDiff{false}; amrex::Real NumDiffCoeff{0.}; diff --git a/Source/ERF.H b/Source/ERF.H index 117aa9b80..e7cedc645 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -336,6 +336,11 @@ public: MultiBlockContainer *m_mbc = nullptr; amrex::Vector > vars_new; amrex::Vector > vars_old; + + // Velocity time averaged field + amrex::Vector> vel_t_avg; + amrex::Vector t_avg_cnt; + #endif std::string pp_prefix {"erf"}; @@ -578,6 +583,10 @@ private: #ifndef ERF_USE_MULTIBLOCK amrex::Vector > vars_new; amrex::Vector > vars_old; + + // Velocity time averaged field + amrex::Vector> vel_t_avg; + amrex::Vector t_avg_cnt; #endif amrex::Vector > > > mri_integrator_mem; @@ -771,6 +780,8 @@ private: "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb", "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac", + // Time averaged velocity + "u_t_avg", "v_t_avg", "w_t_avg", "umag_t_avg", // eddy viscosity "Kmv","Kmh", // eddy diffusivity of heat diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 51ab507bc..e31ec67df 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -263,6 +263,10 @@ ERF::ERF () // Qv prim for MOST Qv_prim.resize(nlevs_max); + // Time averaged velocity field + vel_t_avg.resize(nlevs_max); + t_avg_cnt.resize(nlevs_max); + // Initialize tagging criteria for mesh refinement refinement_criteria_setup(); @@ -321,7 +325,7 @@ ERF::Evolve () cur_time += dt[0]; Print() << "Coarse STEP " << step+1 << " ends." << " TIME = " << cur_time - << " DT = " << dt[0] << std::endl; + << " DT = " << dt[0] << std::endl; post_timestep(step, cur_time, dt[0]); @@ -902,6 +906,16 @@ ERF::InitData () for (int lev = 0; lev <= finest_level; ++lev) micro->Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]); } + // Fill time averaged velocities before first plot file + if (solverChoice.time_avg_vel) { + for (int lev = 0; lev <= finest_level; ++lev) { + Time_Avg_Vel_atCC(dt[lev], t_avg_cnt[lev], vel_t_avg[lev].get(), + vars_new[lev][Vars::xvel], + vars_new[lev][Vars::yvel], + vars_new[lev][Vars::zvel]); + } + } + // check for additional plotting variables that are available after particle containers // are setup. const std::string& pv1 = "plot_vars_1"; appendPlotVariables(pv1,plot_var_names_1); @@ -1847,6 +1861,10 @@ ERF::ERF (const RealBox& rb, int max_level_in, // Theta prim for MOST Theta_prim.resize(nlevs_max); + // Time averaged velocity field + vel_t_avg.resize(nlevs_max); + t_avg_cnt.resize(nlevs_max); + // Initialize tagging criteria for mesh refinement refinement_criteria_setup(); diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 316217410..f0cd5cf53 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -158,6 +158,21 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, rV_new[lev].setVal(3.4e22); rW_new[lev].setVal(5.6e23); + // ******************************************************************************************** + // These are just time averaged fields for diagnostics + // ******************************************************************************************** + + // NOTE: We are not completing a fillpach call on the time averaged data; + // which would copy on intersection and interpolate from coarse. + // Therefore, we are restarting the averaging when the ba changes, + // this may give poor statistics for dynamic mesh refinment. + vel_t_avg[lev] = nullptr; + if (solverChoice.time_avg_vel) { + vel_t_avg[lev] = std::make_unique(ba, dm, 4, 0); // Each vel comp and the mag + vel_t_avg[lev]->setVal(0.0); + t_avg_cnt[lev] = 0.0; + } + // ******************************************************************************************** // Initialize flux registers whenever we create/re-create a level // ******************************************************************************************** @@ -169,7 +184,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, advflux_reg[lev] = new YAFluxRegister(ba , grids[lev-1], dm , dmap[lev-1], geom[lev], geom[lev-1], - ref_ratio[lev-1], lev, ncomp_reflux); + ref_ratio[lev-1], lev, ncomp_reflux); } } @@ -177,15 +192,15 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // Define Theta_prim storage if using MOST BC // ******************************************************************************************** if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { - Theta_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); - if (solverChoice.moisture_type != MoistureType::None) { - Qv_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); - } else { - Qv_prim[lev] = nullptr; - } + Theta_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); + if (solverChoice.moisture_type != MoistureType::None) { + Qv_prim[lev] = std::make_unique(ba,dm,1,IntVect(ngrow_state,ngrow_state,0)); + } else { + Qv_prim[lev] = nullptr; + } } else { - Theta_prim[lev] = nullptr; - Qv_prim[lev] = nullptr; + Theta_prim[lev] = nullptr; + Qv_prim[lev] = nullptr; } // ******************************************************************************************** diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 23be8b4fa..43f8c7dff 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -381,7 +381,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp //******************************************************************************************** // This allocates all kinds of things, including but not limited to: solution arrays, - // terrain arrays and metrics,a nd base state. + // terrain arrays and metrics, and base state. // ******************************************************************************************* init_stuff(lev, ba, dm, temp_lev_new, temp_lev_old); diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 253d353ce..c78dab745 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -749,6 +749,80 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } + if (solverChoice.time_avg_vel) { + if (containerHasElement(plot_var_names, "u_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,0) / norm; + }); + } + mf_comp ++; + } + + if (containerHasElement(plot_var_names, "v_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,1) / norm; + }); + } + mf_comp ++; + } + + if (containerHasElement(plot_var_names, "w_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,2) / norm; + }); + } + mf_comp ++; + } + + if (containerHasElement(plot_var_names, "umag_t_avg")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& data = vel_t_avg[lev]->array(mfi); + const Real norm = t_avg_cnt[lev]; + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + derdat(i ,j ,k, mf_comp) = data(i,j,k,3) / norm; + }); + } + mf_comp ++; + } + } + if (containerHasElement(plot_var_names, "Kmv")) { MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Mom_v,mf_comp,1,0); mf_comp ++; diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index d8b901e03..cc2f4694d 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -179,4 +179,11 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) {time, time + dt_lev}); } } + + // *********************************************************************************************** + // 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); + } } diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index f0b87ba00..2100badf8 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -1,7 +1,3 @@ -CEXE_sources += MomentumToVelocity.cpp -CEXE_sources += VelocityToMomentum.cpp -CEXE_sources += InteriorGhostCells.cpp - CEXE_headers += TerrainMetrics.H CEXE_headers += Microphysics_Utils.H CEXE_headers += TileNoZ.H @@ -11,7 +7,6 @@ CEXE_headers += Interpolation_WENO.H CEXE_headers += Interpolation_WENO_Z.H CEXE_headers += Interpolation.H CEXE_headers += Interpolation_1D.H -CEXE_sources += TerrainMetrics.cpp CEXE_headers += ParFunctions.H @@ -19,6 +14,12 @@ CEXE_headers += Sat_methods.H CEXE_headers += Water_vapor_saturation.H CEXE_headers += DirectionSelector.H +CEXE_sources += MomentumToVelocity.cpp +CEXE_sources += VelocityToMomentum.cpp +CEXE_sources += InteriorGhostCells.cpp +CEXE_sources += TerrainMetrics.cpp +CEXE_sources += Time_Avg_Vel.cpp + ifeq ($(USE_POISSON_SOLVE),TRUE) CEXE_sources += ERF_PoissonSolve.cpp CEXE_sources += ERF_PoissonSolve_tb.cpp diff --git a/Source/Utils/Time_Avg_Vel.cpp b/Source/Utils/Time_Avg_Vel.cpp new file mode 100644 index 000000000..5121ede60 --- /dev/null +++ b/Source/Utils/Time_Avg_Vel.cpp @@ -0,0 +1,47 @@ +#include + +using namespace amrex; + +/* + * Accumulate time averaged velocity fields + */ +void +Time_Avg_Vel_atCC (const Real& dt, + Real& t_avg_cnt, + MultiFab* vel_t_avg, + MultiFab& xvel, + MultiFab& yvel, + MultiFab& zvel) +{ + // Augment the counter + t_avg_cnt += dt; + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*(vel_t_avg),TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // CC tilebox + Box tbx = mfi.tilebox(); + + // Velocity on faces + const Array4& velx = xvel.array(mfi); + const Array4& vely = yvel.array(mfi); + const Array4& velz = zvel.array(mfi); + + // Time average at CC + Array4 vel_t_avg_arr = vel_t_avg->array(mfi); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real u_cc = 0.5 * ( velx(i,j,k) + velx(i+1,j ,k ) ); + Real v_cc = 0.5 * ( vely(i,j,k) + vely(i ,j+1,k ) ); + Real w_cc = 0.5 * ( velz(i,j,k) + velz(i ,j ,k+1) ); + Real umag_cc = std::sqrt(u_cc*u_cc + v_cc*v_cc + w_cc*w_cc); + vel_t_avg_arr(i,j,k,0) += u_cc; + vel_t_avg_arr(i,j,k,1) += v_cc; + vel_t_avg_arr(i,j,k,2) += w_cc; + vel_t_avg_arr(i,j,k,3) += umag_cc; + }); + } +} diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index 4f73531a1..42c25ad1b 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -109,6 +109,17 @@ fine_compute_interior_ghost_rhs (const amrex::Real& time, amrex::Vector& S_rhs_f, amrex::Vector& S_data_f); +/* + * Accumulate time averaged velocity fields + */ +void +Time_Avg_Vel_atCC (const amrex::Real& dt, + amrex::Real& t_avg_cnt, + amrex::MultiFab* vel_t_avg, + amrex::MultiFab& xvel, + amrex::MultiFab& yvel, + amrex::MultiFab& zvel); + /** * Zero RHS in the set region *