ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_MOSTAverage.H
Go to the documentation of this file.
1 #ifndef ERF_MOSTAverage_H
2 #define ERF_MOSTAverage_H
3 
4 #include <AMReX_Gpu.H>
5 #include <AMReX_FArrayBox.H>
6 #include <AMReX_MultiFab.H>
7 #include <AMReX_iMultiFab.H>
8 #include <AMReX_ParmParse.H>
9 #include <ERF_IndexDefines.H>
10 #include <ERF_TerrainMetrics.H>
11 
12 class MOSTAverage {
13 public:
14  explicit MOSTAverage (amrex::Vector<amrex::Geometry> geom,
15  amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
16  amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
17  amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim,
18  amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qr_prim,
19  amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd);
20 
21  // MOSTAverage() = default;
23  {}
24 
25  // Declare a default move constructor so we ensure the destructor is
26  // not called when we return an object of this class by value
27  MOSTAverage (MOSTAverage&&) noexcept = default;
28 
29  // Delete the move assignment operator
30  MOSTAverage& operator=(MOSTAverage&& other) noexcept = delete;
31 
32  // Delete the copy constructor
33  MOSTAverage (const MOSTAverage& other) = delete;
34 
35  // Delete the copy assignment operator
36  MOSTAverage& operator=(const MOSTAverage& other) = delete;
37 
38  // Reset the pointers to field MFs
39  void update_field_ptrs (int lev,
40  amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
41  amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
42  amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim,
43  amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qr_prim);
44 
45  // Update the rotated fields
46  void set_rotated_fields (int lev);
47 
48  // Compute ncells per plane
50 
51  // Compute ncells per plane
53  {m_ncell_region = (2 * m_radius + 1) * (2 * m_radius + 1) * (2 * m_radius + 1);}
54 
55  // Populate a 2D iMF k_indx (w/o terrain)
56  void set_k_indices_N ();
57 
58  // Populate a 2D iMF k_indx (w/ terrain)
59  void set_k_indices_T ();
60 
61  // Populate all 2D iMFs ijk_indx (w/ terrain)
62  void set_norm_indices_T ();
63 
64  // Populate positions (w/ terrain & norm vector & interpolation)
65  void set_z_positions_T ();
66 
67  // Populate positions (w/ terrain & norm vector & interpolation)
68  void set_norm_positions_T ();
69 
70  // Driver for the different average policies
71  void compute_averages (int lev);
72 
73  // Fill averages for policy::plane
74  void compute_plane_averages (int lev);
75 
76  // Fill averages for policy::point
77  void compute_region_averages (int lev);
78 
79  // Write k indices
80  void write_k_indices (int lev);
81 
82  // Write ijk indices
83  void write_norm_indices (int lev);
84 
85  // Write XZ planar positions
86  void write_xz_positions (int lev, int j);
87 
88  // Write averages on 2D mf
89  void write_averages (int lev);
90 
91  // Get pointer to the 2D mf of averages
92  [[nodiscard]] const amrex::MultiFab* get_average (int lev, int comp) const { return m_averages[lev][comp].get(); }
93 
94  // Get z_ref (may be computed from specified k_indx)
95  [[nodiscard]] amrex::Real get_zref () const { return m_zref; }
96 
97  /**
98  * Function to compute trilinear interpolation with terrain.
99  *
100  * @param[in] xp X-position
101  * @param[in] yp Y-position
102  * @param[out] interp_vals Values interpolated
103  * @param[in] interp_array Array to interpolate on
104  * @param[in] z_arr Physical heights
105  * @param[in] plo Problem lower bounds
106  * @param[in] dxi Inverse cell size array
107  * @param[in] interp_comp Number of components to interpolate
108  */
109  AMREX_GPU_HOST_DEVICE AMREX_INLINE
110  static void trilinear_interp_T (const amrex::Real& xp,
111  const amrex::Real& yp,
112  const amrex::Real& zp,
113  amrex::Real* interp_vals,
114  amrex::Array4<amrex::Real const> const& interp_array,
115  amrex::Array4<amrex::Real const> const& z_arr,
116  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& plo,
117  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dxi,
118  const int interp_comp)
119  {
120  // Search to get z/k
121  bool found = false;
122  amrex::Real zval= 0.0;
123  int kmax = ubound(z_arr).z;
124  int i_new = (int) (xp * dxi[0] - 0.5);
125  int j_new = (int) (yp * dxi[1] - 0.5);
126  amrex::Real z_target = zp;
127  for (int lk(0); lk<kmax; ++lk) {
128  amrex::Real z_lo = 0.25 * ( z_arr(i_new,j_new ,lk ) + z_arr(i_new+1,j_new ,lk )
129  + z_arr(i_new,j_new+1,lk ) + z_arr(i_new+1,j_new+1,lk ) );
130  amrex::Real z_hi = 0.25 * ( z_arr(i_new,j_new ,lk+1) + z_arr(i_new+1,j_new ,lk+1)
131  + z_arr(i_new,j_new+1,lk+1) + z_arr(i_new+1,j_new+1,lk+1) );
132  if (z_target > z_lo && z_target < z_hi){
133  found = true;
134  zval = (amrex::Real) lk + ((z_target - z_lo) / (z_hi - z_lo)) + 0.5;
135  break;
136  }
137  }
138 
139  amrex::ignore_unused(found);
140  AMREX_ASSERT_WITH_MESSAGE(found, "MOSTAverage: Height above terrain not found, try increasing z_ref!");
141 
142  const amrex::RealVect lx((xp - plo[0])*dxi[0] + 0.5,
143  (yp - plo[1])*dxi[1] + 0.5,
144  zval);
145 
146  const amrex::IntVect ijk = lx.floor();
147 
148  int i = ijk[0]; int j = ijk[1]; int k = ijk[2];
149 
150  // Weights
151  const amrex::RealVect sx_hi = lx - ijk;
152  const amrex::RealVect sx_lo = 1.0 - sx_hi;
153 
154  for (int n = 0; n < interp_comp; n++)
155  interp_vals[n] = sx_lo[0]*sx_lo[1]*sx_lo[2]*interp_array(i-1, j-1, k-1,n) +
156  sx_lo[0]*sx_lo[1]*sx_hi[2]*interp_array(i-1, j-1, k ,n) +
157  sx_lo[0]*sx_hi[1]*sx_lo[2]*interp_array(i-1, j , k-1,n) +
158  sx_lo[0]*sx_hi[1]*sx_hi[2]*interp_array(i-1, j , k ,n) +
159  sx_hi[0]*sx_lo[1]*sx_lo[2]*interp_array(i , j-1, k-1,n) +
160  sx_hi[0]*sx_lo[1]*sx_hi[2]*interp_array(i , j-1, k ,n) +
161  sx_hi[0]*sx_hi[1]*sx_lo[2]*interp_array(i , j , k-1,n) +
162  sx_hi[0]*sx_hi[1]*sx_hi[2]*interp_array(i , j , k ,n);
163  }
164 
165 protected:
166 
167  // Passed through constructor
168  //--------------------------------------------
169  const amrex::Vector<amrex::Geometry> m_geom; // Geometry at each level
170  amrex::Vector<amrex::Vector<amrex::MultiFab*>> m_fields; // Ptr to fields to be averaged
171  amrex::Vector<amrex::MultiFab*> m_z_phys_nd; // Ptr to terrain height coords
172 
173  // General vars for multiple or all policies
174  //--------------------------------------------
175  int m_nvar{6}; // 6 fields for U/V/T/Qv/Qr/W
176  int m_navg{6}; // 6 averages for U/V/T/Qv/Tv/Umag
177  int m_maxlev{0}; // Total number of levels
178  int m_policy{0}; // Policy for type of averaging
179  bool m_rotate{false}; // Do vector rotations for terrain?
180  amrex::Real m_zref{10.0}; // Height above surface for MOST BC
181  std::string m_pp_prefix {"erf"}; // ParmParse prefix
182  amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_x_pos; // Ptr to 2D mf to hold x position (maxlev)
183  amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_y_pos; // Ptr to 2D mf to hold y position (maxlev)
184  amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_z_pos; // Ptr to 2D mf to hold z position (maxlev)
185  amrex::Vector<std::unique_ptr<amrex::iMultiFab>> m_i_indx; // Ptr to 2D imf to hold i indices (maxlev)
186  amrex::Vector<std::unique_ptr<amrex::iMultiFab>> m_j_indx; // Ptr to 2D imf to hold j indices (maxlev)
187  amrex::Vector<std::unique_ptr<amrex::iMultiFab>> m_k_indx; // Ptr to 2D imf to hold k indices (maxlev)
188  amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> m_averages; // Ptr to 2D mf to hold averages (maxlev,navg)
189  amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>> m_rot_fields; // Rotated field data
190 
191  // Vars for planar average policy
192  //--------------------------------------------
193  amrex::Vector<amrex::Vector<int>> m_ncell_plane; // Number of cells in plane (maxlev,navg)
194  amrex::Vector<amrex::Vector<amrex::Real>> m_plane_average; // Plane avgs (maxlev,navg)
195 
196  // Vars for point/region average policy
197  //--------------------------------------------
198  int m_radius{0}; // Radius around k_index
199  int m_ncell_region{1}; // Number of cells in local region
200  amrex::Vector<int> m_k_in; // Specified k_index for region avg (maxlev)
201 
202  // Vars for normal vector policy
203  //--------------------------------------------
204  bool m_interp{false}; // Do interpolation on destination?
205  bool m_norm_vec{false}; // Use normal vector to find IJK?
206 
207  // Time average w/ exponential filter fun
208  //--------------------------------------------
209  bool m_t_avg{false}; // Flag to do moving average in time
210  amrex::Vector<int> m_t_init; // Flag to specify if averages are initialized
211  amrex::Real m_time_window{1.0e-16}; // Width of the exp filter function
212  amrex::Real m_fact_new, m_fact_old; // Time average factors for new and old means
213 
214  // Surface velocity correction
215  //--------------------------------------------
216  bool include_subgrid_vel = false;
217  amrex::Vector<amrex::Real> m_Vsg; // Subgrid velocity scale (Mahrt & Sun 1995 MWR)
218 };
219 #endif
Definition: ERF_MOSTAverage.H:12
int m_navg
Definition: ERF_MOSTAverage.H:176
bool m_t_avg
Definition: ERF_MOSTAverage.H:209
int m_policy
Definition: ERF_MOSTAverage.H:178
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_averages
Definition: ERF_MOSTAverage.H:188
MOSTAverage(amrex::Vector< amrex::Geometry > geom, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_old, 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< std::unique_ptr< amrex::MultiFab >> &z_phys_nd)
void set_k_indices_T()
Definition: ERF_MOSTAverage.cpp:426
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_y_pos
Definition: ERF_MOSTAverage.H:183
void update_field_ptrs(int lev, amrex::Vector< amrex::Vector< amrex::MultiFab >> &vars_old, 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)
Definition: ERF_MOSTAverage.cpp:228
void set_region_normalization()
Definition: ERF_MOSTAverage.H:52
void write_xz_positions(int lev, int j)
Definition: ERF_MOSTAverage.cpp:1401
void set_norm_indices_T()
Definition: ERF_MOSTAverage.cpp:483
int m_radius
Definition: ERF_MOSTAverage.H:198
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_i_indx
Definition: ERF_MOSTAverage.H:185
void set_rotated_fields(int lev)
Definition: ERF_MOSTAverage.cpp:247
amrex::Vector< amrex::MultiFab * > m_z_phys_nd
Definition: ERF_MOSTAverage.H:171
amrex::Vector< int > m_t_init
Definition: ERF_MOSTAverage.H:210
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_x_pos
Definition: ERF_MOSTAverage.H:182
void set_plane_normalization()
Definition: ERF_MOSTAverage.cpp:306
amrex::Real get_zref() const
Definition: ERF_MOSTAverage.H:95
amrex::Vector< amrex::Real > m_Vsg
Definition: ERF_MOSTAverage.H:217
void write_k_indices(int lev)
Definition: ERF_MOSTAverage.cpp:1313
void set_norm_positions_T()
Definition: ERF_MOSTAverage.cpp:602
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > m_rot_fields
Definition: ERF_MOSTAverage.H:189
void compute_averages(int lev)
Definition: ERF_MOSTAverage.cpp:664
bool m_rotate
Definition: ERF_MOSTAverage.H:179
amrex::Vector< amrex::Vector< amrex::Real > > m_plane_average
Definition: ERF_MOSTAverage.H:194
amrex::Vector< std::unique_ptr< amrex::MultiFab > > m_z_pos
Definition: ERF_MOSTAverage.H:184
void set_k_indices_N()
Definition: ERF_MOSTAverage.cpp:367
amrex::Vector< amrex::Vector< int > > m_ncell_plane
Definition: ERF_MOSTAverage.H:193
std::string m_pp_prefix
Definition: ERF_MOSTAverage.H:181
bool m_norm_vec
Definition: ERF_MOSTAverage.H:205
int m_nvar
Definition: ERF_MOSTAverage.H:175
amrex::Real m_fact_new
Definition: ERF_MOSTAverage.H:212
void set_z_positions_T()
Definition: ERF_MOSTAverage.cpp:556
const amrex::MultiFab * get_average(int lev, int comp) const
Definition: ERF_MOSTAverage.H:92
void compute_region_averages(int lev)
Definition: ERF_MOSTAverage.cpp:975
int m_ncell_region
Definition: ERF_MOSTAverage.H:199
void compute_plane_averages(int lev)
Definition: ERF_MOSTAverage.cpp:690
void write_norm_indices(int lev)
Definition: ERF_MOSTAverage.cpp:1352
int m_maxlev
Definition: ERF_MOSTAverage.H:177
void write_averages(int lev)
Definition: ERF_MOSTAverage.cpp:1431
bool include_subgrid_vel
Definition: ERF_MOSTAverage.H:216
~MOSTAverage()
Definition: ERF_MOSTAverage.H:22
amrex::Real m_time_window
Definition: ERF_MOSTAverage.H:211
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_j_indx
Definition: ERF_MOSTAverage.H:186
MOSTAverage(MOSTAverage &&) noexcept=default
amrex::Vector< amrex::Vector< amrex::MultiFab * > > m_fields
Definition: ERF_MOSTAverage.H:170
amrex::Vector< int > m_k_in
Definition: ERF_MOSTAverage.H:200
amrex::Real m_fact_old
Definition: ERF_MOSTAverage.H:212
bool m_interp
Definition: ERF_MOSTAverage.H:204
const amrex::Vector< amrex::Geometry > m_geom
Definition: ERF_MOSTAverage.H:169
AMREX_GPU_HOST_DEVICE static AMREX_INLINE void trilinear_interp_T(const amrex::Real &xp, const amrex::Real &yp, const amrex::Real &zp, amrex::Real *interp_vals, amrex::Array4< amrex::Real const > const &interp_array, amrex::Array4< amrex::Real const > const &z_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &plo, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxi, const int interp_comp)
Definition: ERF_MOSTAverage.H:110
amrex::Real m_zref
Definition: ERF_MOSTAverage.H:180
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > m_k_indx
Definition: ERF_MOSTAverage.H:187
Definition: ERF_ConsoleIO.cpp:12