Time avg velocity output (#1582)

* Add time average of velocity fields for output; supports industry stakeholders.

* Update the docs.

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

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 |
 |                             |                  |
 |                             |                  |

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

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.};

Source/ERF.H

@@ -336,6 +336,11 @@ public:
     MultiBlockContainer *m_mbc = nullptr;
     amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_new;
     amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_old;
+
+    // Velocity time averaged field
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> vel_t_avg;
+    amrex::Vector<int> t_avg_cnt;
+
 #endif
     std::string pp_prefix {"erf"};
 
@@ -578,6 +583,10 @@ private:
 #ifndef ERF_USE_MULTIBLOCK
     amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_new;
     amrex::Vector<amrex::Vector<amrex::MultiFab> > vars_old;
+
+    // Velocity time averaged field
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> vel_t_avg;
+    amrex::Vector<amrex::Real> t_avg_cnt;
 #endif
     amrex::Vector<std::unique_ptr<MRISplitIntegrator<amrex::Vector<amrex::MultiFab> > > > 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

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();
 

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<MultiFab>(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,23 +184,23 @@ 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);
         }
     }
 
     // ********************************************************************************************
     // 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<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
-      if (solverChoice.moisture_type != MoistureType::None) {
-          Qv_prim[lev]    = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
-      } else {
-          Qv_prim[lev]    = nullptr;
-      }
+        Theta_prim[lev] = std::make_unique<MultiFab>(ba,dm,1,IntVect(ngrow_state,ngrow_state,0));
+        if (solverChoice.moisture_type != MoistureType::None) {
+            Qv_prim[lev]    = std::make_unique<MultiFab>(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;
     }
 
     // ********************************************************************************************

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);
 

Source/IO/Plotfile.cpp

@@ -749,6 +749,80 @@ ERF::WritePlotFile (int which, Vector<std::string> 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<Real>& derdat = mf[lev].array(mfi);
+                    const Array4<Real>& 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<Real>& derdat = mf[lev].array(mfi);
+                    const Array4<Real>& 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<Real>& derdat = mf[lev].array(mfi);
+                    const Array4<Real>& 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<Real>& derdat = mf[lev].array(mfi);
+                    const Array4<Real>& 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 ++;

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);
+    }
 }