Box perturbation multi-face/corner generalization (erf-model#1716)

* fix merge conflict, OPENMP compile time Error, and trailing white space * git push --set-upstream origin July23_perturbation_box_method * file rename. Moving to perlmutter to do GPU test + validations * added in generalization for source/direct perturbation using Richardson number scaling. Added net-zero energy check for box perturbation method. modified some tags in inputs file. Made new functions in TurbPertStruct.H to clean up code. Continue working on validation, multiface initialization, update .rst file * fixed codespell and trailing white space * added conditions to trigger scale ignoring perturbation through gravity term * fixed CI workflow error. * fixed last CI error * fix last CI bug of implicit capture * added multi-face multi-corner generalization to inflow perturbation method. * code spell fix * all input file change for number of integer entry to perturbation_direction * added multi-face multi-corner generalization to inflow perturbation method + updated documentation * added multi-face multi-corner generalization to inflow perturbation method + updated documentation * small modifications to input files. begin testing * moved amp_calc into private member. Now TurbPertStruct.H is fully cleaned up for review. * fixed trailing whitespace * fixed host/device logic issue, and normalizing by 0 issue on host * mergeing with develope * correcting Perlmutter CUDA compile error * merge once more + GNUmakefile edits * got compressible_direct_intputs to give good magnitudes --------- Co-authored-by: Ting-Hsuan (Dustin) Ma <[email protected]> Co-authored-by: Aaron M. Lattanzi <[email protected]>
AMLattanzi · Aug 1, 2024 · fa6acad · fa6acad
1 parent fd1673c
commit fa6acad
Show file tree

Hide file tree

Showing 10 changed files with 498 additions and 261 deletions.
diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst
@@ -242,7 +242,7 @@ The sponge data is input as a text file with 3 columns containing :math:`z, u, v
 Inflow turbulence generation
 ---------------------------
 
-ERF provides the capability to apply a perturbation zone at the inflow domain boundary to mechanically trip turbulence into the domain.
+ERF provides the capability to apply a perturbation zone at the inflow domain boundary to mechanically trip turbulence within the domain. The current version of the turbulence generation techniques allows for `x` and `y` (horizontal) direction perturbations.
 
 .. |PBinflw| image:: figures/PBIllustration.png
            :width: 600
@@ -257,109 +257,99 @@ ERF provides the capability to apply a perturbation zone at the inflow domain bo
    |  Image taken from `DeLeon et al. (2018)`            |
    +-----------------------------------------------------+
 
-Two different types of perturbation are currently available, ``source``, adopted from `DeLeon et al. (2018)`_
+Two different types of perturbation methods are currently available, ``source`` and ``direct``. Both methods uses the formulation introduced by `DeLeon et al. (2018)`_. The ``source`` option applies the perturbation amplitude range, :math:`\pm \Phi_{PB}`, to each cell within the perturbation box as a source term. Conversely, the ``direct`` option applies the calculated temperature difference directly onto the :math:`\rho \theta` field.
 
 ..
-   _`DeLeon et al. (2018)`: https://doi.org/10.2514/1.J057245
+  _`DeLeon et al. (2018)`: https://doi.org/10.2514/1.J057245
 
-and ``direct``, adopted from `Munoz-Esparza et al. (2015)`_. The ``source`` option applies the perturbation amplitude range, :math:`\pm \Phi_{PB}`, to each cell within the perturbation box as a source term. It's important to note that while this perturbation starts as white noise, it becomes colored noise due to the eddy viscosity turbulence closure. Conversely, the ``direct`` option applies the calculated temperature difference directly onto the :math:`\rho \theta field`.
+The perturbation update interval of the individual perturbation box is determined by,
 
-The current implementation only supports West and South face perturbations, specified by ``erf.perturbation_direction``, where the 3 integer inputs represent the `x`, `y`, and `z` directions, respectively. The flow perturbation method requires the dimensions of an individual box input through ``erf.perturbation_box_dim``, with 3 integer inputs representing :math:`nx_{pb}`, :math:`ny_{pb}`, and :math:`nz_{pb}`, respectively. Following the guidance of `Ma and Senocak (2023)`_,
+.. math::
+
+   \frac{t_p \langle U(z) \rangle_{PB}}{D_{PB}} = 1,
+
+The change in the perturbation amplitude is defined as,
+
+.. math::
+
+   {Ri}_{PB} = \frac{g \beta \Delta \overline{\phi} H_{PB}}{{\langle U(z) \rangle}^2_{PB}}.
+
+The current implementation only supports West, South, East, and North face perturbations, specified by ``erf.perturbation_direction``, where the six integer inputs represent the west, south, bottom, east, north, and top domain face perturbations, respectively. Note that while the top and bottom options are included, triggering these will abort the program. In addition to the direction of perturbation, the flow perturbation method requires the dimensions of an individual box, specified through erf.perturbation_box_dim, with three integer inputs representing :math:`{Nx}_{pb}`, :math:`{Ny}_{pb}`, and :math:`{Nz}_{pb}`, respectively. Following the guidance of `Ma and Senocak (2023)`_, the general rule of thumb is to use :math:`H_{PB} = 1/8 \delta` as the height of the perturbation box, where :math:`\delta` is the boundary layer height. The length of the box in the x-direction should be :math:`L_{PB} = 2H_{PB}`. Depending on the direction of the bulk flow, the width of the box in the y-direction should be defined as :math:`W_{PB} = L_{PB} \tan{(\theta_{inflow})}`.
 
 .. _`Ma and Senocak (2023)`: https://link.springer.com/article/10.1007/s10546-023-00786-1
 
-the general rule of thumb is to use :math:`H_{PB} = 1/8 \delta` as the height of the perturbation box, where :math:`\delta` is the boundary layer height. The length of the box in the x-direction should be :math:`L_{PB} = 2H_{PB}`. Depending on the direction of the bulk flow, the width of the box in the y-direction should be defined as :math:`W_{PB} = L_{PB} \tan{\theta_{inflow}}`. Note that the current implementation only accepts ``int`` entries. Therefore, considering the domain size and mesh resolution, the dimensions of a singular box can be determined.
+Note that the current implementation only accepts integer entries. Therefore, considering the domain size and mesh resolution, the dimensions of a single box can be determined.
 
-The specification of the number of layers and the offset into the domain of the perturbation boxes can be made through ``erf.perturbation_layers`` and ``erf.perturbation_offset``, respectively.
+The specification of the number of layers and the offset into the domain of the perturbation boxes can be made through ``erf.perturbation_layers`` and ``erf.perturbation_offset``, respectively. Below is an example of the required input tags to setup a simulation with an inflow perturbation simulation.
 
 ::
 
           erf.inlet_perturbation_type = "source"
 
-          erf.perturbation_direction = 1 0 0
+          erf.perturbation_direction = 1 0 0 0 0 0
           erf.perturbation_box_dims = 8 8 4
           erf.perturbation_layers = 3
-          erf.perturbation_offset = 1
+          erf.perturbation_offset = 5
 
           erf.perturbation_nondimensional = 0.042
           erf.perturbation_T_infinity = 300.0
-          erf.perturbation_T_intensity = 0.1
+          #erf.perturbation_T_intensity = 0.1
+
+The ``erf.perturbation_T_intensity`` tag can be turned on or off by providing a value or commenting it out. When a value is provided, a pseudo-gravity value is used (solved from the Richardson formulation) to normalize the scales of the problem. While this generates quick turbulence, it should be used as a sanity check rather than a runtime strategy for turbulence generation. Additionally, a net-zero energy enforcement is applied over the perturbation boxes to ensure that the synthetic method does not introduce excess energy into the system at each iteration. Below, we provide a detailed description of the two different types of perturbation methods currently existing within ERF.
 
-Before delving into the details, it's important to note that the two methods are interchangeable. While we adhere to the guidelines from the referenced publications, the use of ``direct`` type forcing is not restricted to having unity cell height, nor is ``source`` type forcing limited to boxes. We have generalized the perturbation methods to allow for flexibility in mixing and matching different types of turbulence generation.
+Examples are provided within ``Exec/DevTests/ABL_perturbation_inflow/`` to set up a turbulent open channel flow using inflow/outflow boundary conditions with the aforementioned turbulent inflow generation technique to trigger turbulence within the simulation domain.
 
 Source type forcing
 -------------------
-
-The perturbation update interval is determined by the equation,
-
-.. math::
-
-   \frac{t_p \langle U(z) \rangle_{PB}}{D_{PB}} = 1,
-
-The change in the perturbation is defined as,
-
-.. math::
-
-   {Ri}_{PB} = \frac{g \beta \Delta \overline{\phi} H_{PB}}{{\langle U(z) \rangle}^2_{PB}}.
-
-The magnitude of the perturbation, ignoring the advection and diffusion effects in the transport equation can be made through a proportionality ratio between the update time and change in the box temperature,
+By ignoring the advection and diffusion effects in the transport equation, the amplitude of the peturbation can be made through a proportionality ratio between the update time and change in the box temperature,
 
 .. math::
 
    \Phi_{PB} \propto \frac{\Delta \overline{\phi}}{t_p}
 
-
-and the perturbation amplitude is determined by the equation,
+and the perturbation amplitude is then computed through,
 
 .. math::
 
    \Phi_{PB} = \frac{Ri_{PB} {\langle U(z) \rangle}^3_{PB}}{g \beta D_{PB} H_{PB}}.
 
-The ``source`` type forcing can adopt the box perturbation method by having the following inputs list.
+``source`` type forcing can adopt the box perturbation method by having the following inputs list.
 
 ::
 
           erf.inlet_perturbation_type = "source"
 
-          erf.perturbation_direction = 1 0 0
+          erf.perturbation_direction = 1 0 0 0 0 0
           erf.perturbation_box_dims = 8 8 4
           erf.perturbation_layers = 3
-          erf.perturbation_offset = 1
+          erf.perturbation_offset = 5
 
           erf.perturbation_nondimensional = 0.042 # Ri
           erf.perturbation_T_infinity = 300.0
-          erf.perturbation_T_intensity = 0.1
+          #erf.perturbation_T_intensity = 0.1
 
-The primary purpose of the box perturbation method is not merely to perturb the temperature field in a box format. The key difference between the box and cell perturbation methods lies in how the perturbation is introduced into the temperature source term. Both methods use white noise to generate perturbations, but through the eddy diffusivity scheme and temperature transport, this white noise transforms into colored noise. Although the exact nature of the colored noise is unknown, this method effectively initiates perturbations that develop downstream. While colored noise can be directly added in the cell perturbation method, introducing it directly into the temperature field can cause simulation instability.
+The box perturbation method (BPM) perturbs the temperature field :math:`\rho \theta` in a volume (box) format. Each box computes a perturbation update time and amplitude, then independently update at its respective update interval during runtime. A single perturbation amplitude is seen by the computational cells that falls within this box. A pseudo-random perturbations (ie. white noise) is applied over the range :math: `[-\phi, +\phi]` then introduced to the :math:`\rho \theta` field via source term (:math:`F_{\rho \theta}`). As temperature is transported and through the action of the subgrid-scale (SGS) filter for eddy viscosity, white noise temperature perturbations become colored noise in the velocity field.
+
+Using the source term to perturb the momentum field through buoyancy coupling requires more time compared to the "direct" perturbation method. Turbulence onset can be triggered by adjusting the size of the perturbation box, as the amplitude is heavily influenced by having the two dimensional scales of the perturbation box in the denominator. While there is a general rule of thumb for this method, a sensitivity test should be performed to correctly perturb the momentum field to match the turbulence for a specific problem.
 
 Direct type forcing
 -------------------
+The ``direct`` method can also be used to effectively trip turbulence into the domain. Minute differences exists between the ``direct`` and ``source`` method, with the primary one being the perturbation amplitude is directly solved from the Richardson number relationship as,
 
-The perturbation update interval is determined by the equation,
 .. math::
 
-   \frac{t_p U_{1}}{d_{c}}
+   \Phi_{PB} = \frac{Ri_{PB} {\langle U(z) \rangle}^2_{PB}}{g \beta H_{PB}}.
 
-and the perturbation amplitude is determined by the equation,
-.. math::
-
-   Ec = \frac{{U_g}^2}{\rho c_p \theta_{pm}}.
-
-The ``direct`` type forcing can adopt the cell perturbation method by having the following inputs list.
+To activate the ``direct`` type forcing, set the following inputs list.
 ::
 
           erf.inlet_perturbation_type = "direct"
 
-          erf.perturbation_direction = 1 0 0
-          erf.perturbation_box_dims = 8 8 1
+          erf.perturbation_direction = 1 0 0 0 0 0
+          erf.perturbation_box_dims = 16 16 8
           erf.perturbation_layers = 3
-          erf.perturbation_offset = 1
-
-          erf.perturbation_nondimensional = 0.16 #Ec
-          erf.perturbation_rho_0 = 1.0
-          erf.perturbation_cp = 1250
-
-From `Munoz-Esparza et al. (2015)`_ the choice of the Eckert number is 0.16.
+          erf.perturbation_offset = 5
 
-..
-   _`Munoz-Esparza et al. (2015)`: https://doi.org/10.1063/1.4913572
+          erf.perturbation_nondimensional = 0.042 # Ri
+          erf.perturbation_T_infinity = 300.0
+          #erf.perturbation_T_intensity = 0.1
diff --git a/Exec/DevTests/ABL_perturbation_inflow/GNUmakefile b/Exec/DevTests/ABL_perturbation_inflow/GNUmakefile
@@ -20,13 +20,12 @@ USE_SYCL = FALSE
 
 # Debugging
 #DEBUG = TRUE
-DEBUG = FALSE
 
 TEST = TRUE
 USE_ASSERTION = TRUE
 
 # Incompressible Solver
-USE_POISSON_SOLVE = TRUE
+#USE_POISSON_SOLVE = TRUE
 
 # GNU Make
 Bpack := ./Make.package

diff --git a/Exec/DevTests/ABL_perturbation_inflow/compressible_direct_inputs b/Exec/DevTests/ABL_perturbation_inflow/compressible_direct_inputs
@@ -0,0 +1,86 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+stop_time = 180.0
+#max_step = 1
+
+amrex.fpe_trap_invalid = 1
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+#Larger problem
+geometry.prob_extent =  40   2.5     10
+amr.n_cell           = 256    16     64
+
+# Quick debug problem
+#geometry.prob_extent =  20     5     10
+#amr.n_cell           =  64    16     32
+
+geometry.is_periodic = 0 1 0
+
+xlo.type = "Inflow"
+xhi.type = "Outflow"
+zhi.type = "SlipWall"
+#zlo.type = "NoSlipWall"
+
+zlo.type = "Most"
+erf.most.flux_type = "custom"
+erf.most.ustar  = 0.0395   # z=10.
+erf.most.tstar  = 0. # theta flux
+erf.most.qstar  = 0. # qv    flux
+
+xlo.density =   1.0
+xlo.theta   = 300.0
+#xlo.velocity = 1.0 0.0 0.0
+xlo.dirichlet_file = "input_ReTau395Ana_inflow.txt"
+
+# TIME STEP CONTROL
+erf.fixed_dt = 0.0002
+erf.fixed_mri_dt_ratio = 4
+erf.dynamicViscosity = 0.001
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 100    # timesteps between computing mass
+erf.pert_interval  = 100    # timesteps between perturbation output message XXX
+erf.v              = 0       # verbosity in ERF.cpp XXX
+amr.v              = 0       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0      # maximum level number allowed
+
+# PLOTFILES
+erf.plot_file_1     = plt    # prefix of plotfile name
+erf.plot_per_1      = 0.1
+#erf.plot_int_1      = 5
+erf.plot_vars_1     = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta
+
+# CHECKPOINT FILES
+#erf.check_file      = chk    # root name of checkpoint file
+#erf.check_per       = 1.0
+
+# SOLVER CHOICE
+erf.alpha_T = 0.0
+erf.alpha_C = 0.0
+erf.molec_diff_type = "None"
+erf.les_type        = "Smagorinsky"
+erf.Cs              = 0.1
+erf.use_gravity     = true
+erf.buoyancy_type   = 1
+
+# Turbulent inflow generation
+#erf.perturbation_type = "source"
+erf.perturbation_type = "direct"
+erf.perturbation_direction = 1 0 0 0 0 0
+erf.perturbation_layers = 3
+erf.perturbation_offset = 3
+
+erf.perturbation_box_dims = 16 16 8 
+erf.perturbation_nondimensional = 0.001 # Ri
+erf.perturbation_T_infinity = 300.0
+#erf.perturbation_T_intensity = 0.01
+
+# Initial condition for the entire field
+erf.init_type = "input_sounding"
+erf.input_sounding_file = "input_ReTau395Ana_sounding.txt"
+
+# PROBLEM PARAMETERS
+prob.rho_0 = 1.0
+prob.T_0 = 300.0
diff --git a/.../ABL_perturbation_inflow/toc_inout_inputs → ...rbation_inflow/compressible_source_inputs b/.../ABL_perturbation_inflow/toc_inout_inputs → ...rbation_inflow/compressible_source_inputs
@@ -1,5 +1,5 @@
 # ------------------  INPUTS TO MAIN PROGRAM  -------------------
-stop_time = 60.0
+stop_time = 180.0
 #max_step = 5
 
 erf.no_substepping = 1
@@ -8,8 +8,8 @@ fabarray.mfiter_tile_size = 1024 1024 1024
 
 # PROBLEM SIZE & GEOMETRY
 #Larger problem
-geometry.prob_extent =  40     5     10
-amr.n_cell           = 256    32     64
+geometry.prob_extent =  40   2.5     10
+amr.n_cell           = 256    16     64
 
 # Quick debug problem
 #geometry.prob_extent =  20     5     10
@@ -25,7 +25,6 @@ zhi.type = "SlipWall"
 zlo.type = "Most"
 erf.most.flux_type = "custom"
 erf.most.ustar  = 0.0395   # z=10.
-#erf.most.ustar  = 0.395   # z=1.0
 erf.most.tstar  = 0. # theta flux
 erf.most.qstar  = 0. # qv    flux
 
@@ -35,21 +34,23 @@ xlo.theta   = 300.0
 xlo.dirichlet_file = "input_ReTau395Ana_inflow.txt"
 
 # TIME STEP CONTROL
-erf.cfl = 0.5
+#erf_.cfl = 0.8
+erf.fixed_dt = 0.0001
+erf.fixed_mri_dt_ratio = 4
 erf.dynamicViscosity = 0.001
 
 # DIAGNOSTICS & VERBOSITY
-erf.sum_interval   = 1       # timesteps between computing mass
-erf.pert_interval  = 1       # timesteps between perturbation output message XXX
-erf.v              = 0       # verbosity in ERF.cpp XXX
+erf.sum_interval   = 100     # timesteps between computing mass
+erf.pert_interval  = 100     # timesteps between perturbation output message XXX
+erf.v              = 1       # verbosity in ERF.cpp XXX
 amr.v              = 0       # verbosity in Amr.cpp
 
 # REFINEMENT / REGRIDDING
 amr.max_level       = 0      # maximum level number allowed
 
 # PLOTFILES
 erf.plot_file_1     = plt    # prefix of plotfile name
-erf.plot_per_1      = 0.05
+erf.plot_per_1      = 0.1
 #erf.plot_int_1      = 5
 erf.plot_vars_1     = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta
 
@@ -64,18 +65,18 @@ erf.molec_diff_type = "None"
 erf.les_type        = "Smagorinsky"
 erf.Cs              = 0.1
 erf.use_gravity     = true
-erf.buoyancy_type   = 2
+erf.buoyancy_type   = 1
 
 # Turbulent inflow generation
 erf.perturbation_type = "source"
-erf.perturbation_direction = 1 0 0
+erf.perturbation_direction = 1 0 0 0 0 0
 erf.perturbation_layers = 3
 erf.perturbation_offset = 3
 
 erf.perturbation_box_dims = 8 8 4
-erf.perturbation_nondimensional = 0.042 # Ri
+erf.perturbation_nondimensional = 0.1 # Ri
 erf.perturbation_T_infinity = 300.0
-erf.perturbation_T_intensity = 0.1
+#erf.perturbation_T_intensity = 0.1
 
 # Initial condition for the entire field
 #erf.init_type = "uniform"