1# Auxiliary Functionality 2This section documents functionality that is not apart of the core PDE solver, but is used for other miscellaneous tasks, such as statistics collection or in-situ machine learning. 3 4(aux-statistics)= 5## Statistics Collection 6For scale-resolving simulations (such as LES and DNS), statistics for a simulation are more often useful than time-instantaneous snapshots of the simulation itself. 7To make this process more computationally efficient, averaging in the spanwise direction, if physically correct, can help reduce the amount of simulation time needed to get converged statistics. 8 9First, let's more precisely define what we mean by spanwise average. 10Denote $\langle \phi \rangle$ as the Reynolds average of $\phi$, which in this case would be a average over the spanwise direction and time: 11 12$$ 13\langle \phi \rangle(x,y) = \frac{1}{L_z + (T_f - T_0)}\int_0^{L_z} \int_{T_0}^{T_f} \phi(x, y, z, t) \mathrm{d}t \mathrm{d}z 14$$ 15 16where $z$ is the spanwise direction, the domain has size $[0, L_z]$ in the spanwise direction, and $[T_0, T_f]$ is the range of time being averaged over. 17Note that here and in the code, **we assume the spanwise direction to be in the $z$ direction**. 18 19To discuss the details of the implementation we'll first discuss the spanwise integral, then the temporal integral, and lastly the statistics themselves. 20 21### Spanwise Integral 22The function $\langle \phi \rangle (x,y)$ is represented on a 2-D finite element grid, taken from the full domain mesh itself. 23If isoperiodicity is set, the periodic face is extracted as the spanwise statistics mesh. 24Otherwise the negative z face is used. 25We'll refer to this mesh as the *parent grid*, as for every "parent" point in the parent grid, there are many "child" points in the full domain. 26Define a function space on the parent grid as $\mathcal{V}_p^\mathrm{parent} = \{ \bm v(\bm x) \in H^{1}(\Omega_e^\mathrm{parent}) \,|\, \bm v(\bm x_e(\bm X)) \in P_p(\bm{I}), e=1,\ldots,N_e \}$. 27We enforce that the order of the parent FEM space is equal to the full domain's order. 28 29Many statistics are the product of 2 or more solution functions, which results in functions of degree higher than the parent FEM space, $\mathcal{V}_p^\mathrm{parent}$. 30To represent these higher-order functions on the parent FEM space, we perform an $L^2$ projection. 31Define the spanwise averaged function as: 32 33$$ 34\langle \phi \rangle_z(x,y,t) = \frac{1}{L_z} \int_0^{L_z} \phi(x, y, z, t) \mathrm{d}z 35$$ 36 37where the function $\phi$ may be the product of multiple solution functions and $\langle \phi \rangle_z$ denotes the spanwise average. 38The projection of a function $u$ onto the parent FEM space would look like: 39 40$$ 41\bm M u_N = \int_0^{L_x} \int_0^{L_y} u \psi^\mathrm{parent}_N \mathrm{d}y \mathrm{d}x 42$$ 43where $\bm M$ is the mass matrix for $\mathcal{V}_p^\mathrm{parent}$, $u_N$ the coefficients of the projected function, and $\psi^\mathrm{parent}_N$ the basis functions of the parent FEM space. 44Substituting the spanwise average of $\phi$ for $u$, we get: 45 46$$ 47\bm M [\langle \phi \rangle_z]_N = \int_0^{L_x} \int_0^{L_y} \left [\frac{1}{L_z} \int_0^{L_z} \phi(x,y,z,t) \mathrm{d}z \right ] \psi^\mathrm{parent}_N(x,y) \mathrm{d}y \mathrm{d}x 48$$ 49 50The triple integral in the right hand side is just an integral over the full domain 51 52$$ 53\bm M [\langle \phi \rangle_z]_N = \frac{1}{L_z} \int_\Omega \phi(x,y,z,t) \psi^\mathrm{parent}_N(x,y) \mathrm{d}\Omega 54$$ 55 56We need to evaluate $\psi^\mathrm{parent}_N$ at quadrature points in the full domain. 57To do this efficiently, **we assume and exploit the full domain grid to be a tensor product in the spanwise direction**. 58This assumption means quadrature points in the full domain have the same $(x,y)$ coordinate location as quadrature points in the parent domain. 59This also allows the use of the full domain quadrature weights for the triple integral. 60 61### Temporal Integral/Averaging 62To calculate the temporal integral, we do a running average using left-rectangle rule. 63At the beginning of each simulation, the time integral of a statistic is set to 0, $\overline{\phi} = 0$. 64Periodically, the integral is updated using left-rectangle rule: 65 66$$\overline{\phi}_\mathrm{new} = \overline{\phi}_{\mathrm{old}} + \phi(t_\mathrm{new}) \Delta T$$ 67where $\phi(t_\mathrm{new})$ is the statistic at the current time and $\Delta T$ is the time since the last update. 68When stats are written out to file, this running sum is then divided by $T_f - T_0$ to get the time average. 69 70With this method of calculating the running time average, we can plug this into the $L^2$ projection of the spanwise integral: 71 72$$ 73\bm M [\langle \phi \rangle]_N = \frac{1}{L_z + (T_f - T_0)} \int_\Omega \int_{T_0}^{T_f} \phi(x,y,z,t) \psi^\mathrm{parent}_N \mathrm{d}t \mathrm{d}\Omega 74$$ 75where the integral $\int_{T_0}^{T_f} \phi(x,y,z,t) \mathrm{d}t$ is calculated on a running basis. 76 77 78### Running 79As the simulation runs, it takes a running time average of the statistics at the full domain quadrature points. 80This running average is only updated at the interval specified by `-ts_monitor_spanstats_turbulence_collect_interval` as number of timesteps. 81The $L^2$ projection problem is only solved when statistics are written to file, which is controlled by `-ts_monitor_spanstats_turbulence_interval`. 82Note that the averaging is not reset after each file write. 83The average is always over the bounds $[T_0, T_f]$, where $T_f$ in this case would be the time the file was written at and $T_0$ is the solution time at the beginning of the run. 84 85:::{list-table} Spanwise Turbulent Statistics Runtime Options 86:header-rows: 1 87 88* - Option 89 - Description 90 - Default value 91 92* - `-ts_monitor_spanstats_turbulence_collect_interval` 93 - Number of timesteps between statistics collection 94 - `1` 95 96* - `-ts_monitor_spanstats_turbulence` 97 - Sets the `PetscViewer` for the statistics file writing, such as `cgns:output-%d.cgns` (requires PETSc `--download-cgns`). Also turns the statistics collection on. 98 - 99 100* - `-ts_monitor_spanstats_turbulence_interval` 101 - Number of timesteps between statistics file writing (`-1` means only at end of run) 102 - `-1` 103 104* - `-ts_monitor_spanstats_turbulence_cgns_batch_size` 105 - Number of frames written per CGNS file if the CGNS file name includes a format specifier (`%d`). 106 - `20` 107::: 108 109### Turbulent Statistics 110 111The focus here are those statistics that are relevant to turbulent flow. 112The terms collected are listed below, with the mathematical definition on the left and the label (present in CGNS output files) is on the right. 113 114| Math | Label | 115| ----------------- | -------- | 116| $\langle \rho \rangle$ | MeanDensity | 117| $\langle p \rangle$ | MeanPressure | 118| $\langle p^2 \rangle$ | MeanPressureSquared | 119| $\langle p u_i \rangle$ | MeanPressureVelocity[$i$] | 120| $\langle \rho T \rangle$ | MeanDensityTemperature | 121| $\langle \rho T u_i \rangle$ | MeanDensityTemperatureFlux[$i$] | 122| $\langle \rho u_i \rangle$ | MeanMomentum[$i$] | 123| $\langle \rho u_i u_j \rangle$ | MeanMomentumFlux[$ij$] | 124| $\langle u_i \rangle$ | MeanVelocity[$i$] | 125 126where [$i$] are suffixes to the labels. So $\langle \rho u_x u_y \rangle$ would correspond to MeanMomentumFluxXY. 127This naming convention is chosen to align with the CGNS standard naming style. 128 129To get second-order statistics from these terms, simply use the identity: 130 131$$ 132\langle \phi' \theta' \rangle = \langle \phi \theta \rangle - \langle \phi \rangle \langle \theta \rangle 133$$ 134 135(aux-differential-filtering)= 136## Differential Filtering 137 138There is the option to filter the solution field using differential filtering. 139This was first proposed in {cite}`germanoDiffFilterLES1986`, using an inverse Hemholtz operator. 140The strong form of the differential equation is 141 142$$ 143\overline{\phi} - \nabla \cdot (\beta (\bm{D}\bm{\Delta})^2 \nabla \overline{\phi} ) = \phi 144$$ 145 146for $\phi$ the scalar solution field we want to filter, $\overline \phi$ the filtered scalar solution field, $\bm{\Delta} \in \mathbb{R}^{3 \times 3}$ a symmetric positive-definite rank 2 tensor defining the width of the filter, $\bm{D}$ is the filter width scaling tensor (also a rank 2 SPD tensor), and $\beta$ is a kernel scaling factor on the filter tensor. 147This admits the weak form: 148 149$$ 150\int_\Omega \left( v \overline \phi + \beta \nabla v \cdot (\bm{D}\bm{\Delta})^2 \nabla \overline \phi \right) \,d\Omega 151- \cancel{\int_{\partial \Omega} \beta v \nabla \overline \phi \cdot (\bm{D}\bm{\Delta})^2 \bm{\hat{n}} \,d\partial\Omega} = 152\int_\Omega v \phi \, , \; \forall v \in \mathcal{V}_p 153$$ 154 155The boundary integral resulting from integration-by-parts is crossed out, as we assume that $(\bm{D}\bm{\Delta})^2 = \bm{0} \Leftrightarrow \overline \phi = \phi$ at boundaries (this is reasonable at walls, but for convenience elsewhere). 156 157### Filter Width Tensor, Δ 158For homogenous filtering, $\bm{\Delta}$ is defined as the identity matrix. 159 160:::{note} 161It is common to denote a filter width dimensioned relative to the radial distance of the filter kernel. 162Note here we use the filter *diameter* instead, as that feels more natural (albeit mathematically less convenient). 163For example, under this definition a box filter would be defined as: 164 165$$ 166B(\Delta; \bm{r}) = 167\begin{cases} 1681 & \Vert \bm{r} \Vert \leq \Delta/2 \\ 1690 & \Vert \bm{r} \Vert > \Delta/2 170\end{cases} 171$$ 172::: 173 174For inhomogeneous anisotropic filtering, we use the finite element grid itself to define $\bm{\Delta}$. 175This is set via `-diff_filter_grid_based_width`. 176Specifically, we use the filter width tensor defined in {cite}`prakashDDSGSAnisotropic2022`. 177For finite element grids, the filter width tensor is most conveniently defined by $\bm{\Delta} = \bm{g}^{-1/2}$ where $\bm g = \nabla_{\bm x} \bm{X} \cdot \nabla_{\bm x} \bm{X}$ is the metric tensor. 178 179### Filter Width Scaling Tensor, $\bm{D}$ 180The filter width tensor $\bm{\Delta}$, be it defined from grid based sources or just the homogenous filtering, can be scaled anisotropically. 181The coefficients for that anisotropic scaling are given by `-diff_filter_width_scaling`, denoted here by $c_1, c_2, c_3$. 182The definition for $\bm{D}$ then becomes 183 184$$ 185\bm{D} = 186\begin{bmatrix} 187 c_1 & 0 & 0 \\ 188 0 & c_2 & 0 \\ 189 0 & 0 & c_3 \\ 190\end{bmatrix} 191$$ 192 193In the case of $\bm{\Delta}$ being defined as homogenous, $\bm{D}\bm{\Delta}$ means that $\bm{D}$ effectively sets the filter width. 194 195The filtering at the wall may also be damped, to smoothly meet the $\overline \phi = \phi$ boundary condition at the wall. 196The selected damping function for this is the van Driest function {cite}`vandriestWallDamping1956`: 197 198$$ 199\zeta = 1 - \exp\left(-\frac{y^+}{A^+}\right) 200$$ 201 202where $y^+$ is the wall-friction scaled wall-distance ($y^+ = y u_\tau / \nu = y/\delta_\nu$), $A^+$ is some wall-friction scaled scale factor, and $\zeta$ is the damping coefficient. 203For this implementation, we assume that $\delta_\nu$ is constant across the wall and is defined by `-diff_filter_friction_length`. 204$A^+$ is defined by `-diff_filter_damping_constant`. 205 206To apply this scalar damping coefficient to the filter width tensor, we construct the wall-damping tensor from it. 207The construction implemented currently limits damping in the wall parallel directions to be no less than the original filter width defined by $\bm{\Delta}$. 208The wall-normal filter width is allowed to be damped to a zero filter width. 209It is currently assumed that the second component of the filter width tensor is in the wall-normal direction. 210Under these assumptions, $\bm{D}$ then becomes: 211 212$$ 213\bm{D} = 214\begin{bmatrix} 215 \max(1, \zeta c_1) & 0 & 0 \\ 216 0 & \zeta c_2 & 0 \\ 217 0 & 0 & \max(1, \zeta c_3) \\ 218\end{bmatrix} 219$$ 220 221### Filter Kernel Scaling, β 222While we define $\bm{D}\bm{\Delta}$ to be of a certain physical filter width, the actual width of the implied filter kernel is quite larger than "normal" kernels. 223To account for this, we use $\beta$ to scale the filter tensor to the appropriate size, as is done in {cite}`bullExplicitFilteringExact2016`. 224To match the "size" of a normal kernel to our differential kernel, we attempt to have them match second order moments with respect to the prescribed filter width. 225To match the box and Gaussian filters "sizes", we use $\beta = 1/10$ and $\beta = 1/6$, respectively. 226$\beta$ can be set via `-diff_filter_kernel_scaling`. 227 228### Runtime Options 229 230:::{list-table} Differential Filtering Runtime Options 231:header-rows: 1 232 233* - Option 234 - Description 235 - Default value 236 - Unit 237 238* - `-diff_filter_monitor` 239 - Enable differential filter TSMonitor 240 - `false` 241 - boolean 242 243* - `-diff_filter_grid_based_width` 244 - Use filter width based on the grid size 245 - `false` 246 - boolean 247 248* - `-diff_filter_width_scaling` 249 - Anisotropic scaling for filter width in wall-aligned coordinates (snz) 250 - `1,1,1` 251 - `m` 252 253* - `-diff_filter_kernel_scaling` 254 - Scaling to make differential kernel size equivalent to other filter kernels 255 - `0.1` 256 - `m^2` 257 258* - `-diff_filter_wall_damping_function` 259 - Damping function to use at the wall for anisotropic filtering (`none`, `van_driest`) 260 - `none` 261 - string 262 263* - `-diff_filter_wall_damping_constant` 264 - Constant for the wall-damping function. $A^+$ for `van_driest` damping function. 265 - 25 266 - 267 268* - `-diff_filter_friction_length` 269 - Friction length associated with the flow, $\delta_\nu$. Used in wall-damping functions 270 - 0 271 - `m` 272::: 273 274(aux-in-situ-ml)= 275## *In Situ* Machine-Learning Model Training 276Training machine-learning models normally uses *a priori* (already gathered) data stored on disk. 277This is computationally inefficient, particularly as the scale of the problem grows and the data that is saved to disk reduces to a small percentage of the total data generated by a simulation. 278One way of working around this to to train a model on data coming from an ongoing simulation, known as *in situ* (in place) learning. 279 280This is implemented in the code using [SmartSim](https://www.craylabs.org/docs/overview.html). 281Briefly, the fluid simulation will periodically place data for training purposes into a database that a separate process uses to train a model. 282The database used by SmartSim is [Redis](https://redis.com/modules/redis-ai/) and the library to connect to the database is called [SmartRedis](https://www.craylabs.org/docs/smartredis.html). 283More information about how to utilize this code in a SmartSim configuration can be found on [SmartSim's website](https://www.craylabs.org/docs/overview.html). 284 285To use this code in a SmartSim *in situ* setup, first the code must be built with SmartRedis enabled. 286This is done by specifying the installation directory of SmartRedis using the `SMARTREDIS_DIR` environment variable when building: 287 288``` 289make SMARTREDIS_DIR=~/software/smartredis/install 290``` 291 292### SGS Data-Driven Model *In Situ* Training 293Currently the code is only setup to do *in situ* training for the SGS data-driven model. 294Training data is split into the model inputs and outputs. 295The model inputs are calculated as the same model inputs in the SGS Data-Driven model described {ref}`earlier<sgs-dd-model>`. 296The model outputs (or targets in the case of training) are the subgrid stresses. 297Both the inputs and outputs are computed from a filtered velocity field, which is calculated via {ref}`aux-differential-filtering`. 298The settings for the differential filtering used during training are described in {ref}`aux-differential-filtering`. 299The training will create multiple sets of data per each filter width defined in `-sgs_train_filter_widths`. 300Those scalar filter widths correspond to the scaling correspond to $\bm{D} = c \bm{I}$, where $c$ is the scalar filter width. 301 302The SGS *in situ* training can be enabled using the `-sgs_train_enable` flag. 303Data can be processed and placed into the database periodically. 304The interval between is controlled by `-sgs_train_write_data_interval`. 305There's also the choice of whether to add new training data on each database write or to overwrite the old data with new data. 306This is controlled by `-sgs_train_overwrite_data`. 307 308The database may also be located on the same node as a MPI rank (collocated) or located on a separate node (distributed). 309It's necessary to know how many ranks are associated with each collocated database, which is set by `-smartsim_collocated_database_num_ranks`. 310 311### Runtime Options 312:::{list-table} *In Situ* Machine-Learning Training Runtime Options 313:header-rows: 1 314 315* - Option 316 - Description 317 - Default value 318 - Unit 319 320* - `-sgs_train_enable` 321 - Whether to enable *in situ* training of data-driven SGS model. Require building with SmartRedis. 322 - `false` 323 - boolean 324 325* - `-sgs_train_write_data_interval` 326 - Number of timesteps between writing training data into SmartRedis database 327 - `1` 328 - 329 330* - `-sgs_train_overwrite_data` 331 - Whether new training data should overwrite old data on database 332 - `true` 333 - boolean 334 335* - `-sgs_train_filter_widths` 336 - List of scalar values for different filter widths to calculate for training data 337 - 338 - `m` 339 340* - `-smartsim_collocated_num_ranks` 341 - Number of MPI ranks associated with each collocated database (i.e. ranks per node) 342 - `1` 343 - 344::: 345