Turbulent inflow generation logic patch (#1695)

* Generalized size of boxes. Start implementing CPB. Still need user input + .rst documentation

* Added my own folder in Exec, ABL_perturbation_inflow

* Generalized BPM so it can transform into CPM. Added more user input variables. Renamed BPM and CPM tags to source and direct to keep neutrality. Need to add Ec, tp, and cell_data equations for direction add.

* manually got rid of tabs

* added in calc_tpi_meanMag_firstCell() function and direct field add within TimeIntegration/ERF_Advance.cpp for CPM cases. Need to implement options to swap back and forth on use side between BPM and CPM. Switch gear and start runtime optimization.

* amr.v=0 helps with runtime. Added random numbers for PB_intervals. Added outputs for PB_amplitude/PB_updateTime in seperate file when erf.v=1. Need to find a good way to generalize amplitude. Currently breaks code when amplitude is too large. CONTINUE with CPM debug, fast_rhs_fun debug with PB, and documentation with .rst file.

* Fixed CUDA private member issue. CONTINUE with CPM debug, fast_rhs_fun debug with PB, and documentation with .rst file

* Added output format into file + verbose tags. Fixed parallel issue of boxes. Added in perturbation inflow .rst documents. CONTINUE with CPM debug, fast_rhs_fun debug with PB.

* Added output format into file + verbose tags. Fixed parallel issue of boxes. Added in perturbation inflow .rst documents. CONTINUE with CPM debug, fast_rhs_fun debug with PB.

* Get rid of trailing white space. Added output format into file + verbose tags. Fixed parallel issue of boxes. Added in perturbation inflow .rst documents. CONTINUE with CPM debug, fast_rhs_fun debug with PB.

* update DOI links for build-and-deploy fail test.

* 404 Client Error: Not Found for url error fix.

* 404 Client Error: Not Found for url error fix.

* 404 Client Error: Not Found for url error fix.

* Error 403 with Munoz-Esparza 2015 paper. Last attempt

* Added in perturbation inflow .rst documents. CONTINUE with pb tests. Currently blows up after 1st update interval.

* fix for doc push

* manual rerun for DOI link in docs

* link replacement for syntax error

* manual check for munoz link

* what is going on with these links?

* using link I know that works

* testing DeLeon link

* DeLeon: maybe it's a syntax error?

* DeLeon: Trying the other way of linking

* maybe special character problem?

* commented out links to DeLeon et al. (2018), and Munoz-Esparza et al. (2015). DOI links return Erorror 403, actual link returns Error 404.

* testing perturbation inflow method, qualitativly for now

* Added trigger tags for PBM in ERF.cpp. This should fix compile error in workflow.

* Document update

* resolved erroneous error with pseduo gravity (ignores scales with mechancial perturbation)

* resolved erroneous error with pseduo gravity (ignores scales with mechancial perturbation)

* updating PR

* fixing codespell

* corrected math mode in .rst file

* modifying test files

* got rid of some weird merge conflict HEAD that remained

* Should not be using sounding with incompressible flow, as it is a compressible flow algorithm.

* fix to input tag

* reworked apply_tpi() logic so it is now a per cell within PB assignment. Random assignment happens at the calc_tpi_update() level

* white space fix

* fixed unsed variable error for CI tests

* changed input file. modified python file to generate ReTau180, and generate dirichlet input file. Added more description into documentation.

* init perturbation in initialization to see if turbulence is sustained

* Box perturbation method working with incompressible solver.

* fixed codespell

* fixed white space

---------

Co-authored-by: Ting-Hsuan (Dustin) Ma <[email protected]>

Docs/sphinx_doc/BoundaryConditions.rst

@@ -262,13 +262,13 @@ Two different types of perturbation are currently available, ``source``, adopted
 ..
    _`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, `\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 `\rho \theta field`.
+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 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 `nx_{pb}`, `ny_{pb}`, and `nz_{pb}`, respectively. Following the guidance of `Ma and Senocak (2023)`_,
+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)`_,
 
 .. _`Ma and Senocak (2023)`: https://link.springer.com/article/10.1007/s10546-023-00786-1
 
-the general rule of thumb is to use `H_{PB} = 1/8 \delta` as the height of the perturbation box, where `\delta` is the boundary layer height. The length of the box in the x-direction should be `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 `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.
+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.
 
 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.
 
@@ -306,7 +306,8 @@ The magnitude of the perturbation, ignoring the advection and diffusion effects
 
 .. math::
 
-   \Phi_{PB} \propto \frac{\Delta \overliner{\phi}}{t_p}
+   \Phi_{PB} \propto \frac{\Delta \overline{\phi}}{t_p}
+
 
 and the perturbation amplitude is determined by the equation,
 
@@ -329,6 +330,8 @@ The ``source`` type forcing can adopt the box perturbation method by having the
           erf.perturbation_T_infinity = 300.0
           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.
+
 Direct type forcing
 -------------------
 

Docs/sphinx_doc/Inputs.rst

@@ -563,7 +563,7 @@ List of Parameters
 | Parameter                     | Definition       | Acceptable     | Default        |
 |                               |                  | Values         |                |
 +===============================+==================+================+================+
-| **erf.datalog**               | Output           | Up to four     | NONE           |
+| **erf.data_log**              | Output           | Up to four     | NONE           |
 |                               | filename(s)      | strings        |                |
 +-------------------------------+------------------+----------------+----------------+
 | **erf.profile_int**           | Interval (number)| Integer        | -1             |

Exec/DevTests/ABL_perturbation_inflow/GNUmakefile

@@ -19,12 +19,14 @@ USE_HIP  = FALSE
 USE_SYCL = FALSE
 
 # Debugging
-DEBUG = TRUE
+#DEBUG = TRUE
+DEBUG = FALSE
 
 TEST = TRUE
 USE_ASSERTION = TRUE
 
-#USE_POISSON_SOLVE = TRUE
+# Incompressible Solver
+USE_POISSON_SOLVE = TRUE
 
 # GNU Make
 Bpack := ./Make.package

Exec/DevTests/ABL_perturbation_inflow/README

@@ -5,3 +5,9 @@ boundary condition possibly specified by Monin Obukhov Similarity Theory (MOST).
 This version of the ABL problem initializes the data from an input sounding file,
 and initializes turbulence at the inflow. The target is to transition turbulence
 from coarse to fine grid interface.
+
+To run this test, please run the python script to generate the correct data input
+for ERF through, `python3 generateLogLawProfile.py`. Within the python file, please
+edit `Re_tau`, `height_end`, and `kinematic_viscosity` values to the target simulation.
+Once this is done, double check that the input file for ERF corresponds with the 
+files created by the python script.

Exec/DevTests/ABL_perturbation_inflow/generateLogLawProfile.py

@@ -0,0 +1,49 @@
+import numpy as np
+
+# Given parameters
+Re_tau = 395
+kinematic_viscosity = 0.001
+height_start = 1e-2
+height_end = 10
+#height_end = 1.0
+num_points = 256
+kappa = 0.41  # Von Kármán constant
+B = 5.2
+
+# Calculate friction velocity (u_tau) based on Re_tau, kinematic viscosity, and height_end
+u_tau = (Re_tau * kinematic_viscosity) / height_end
+print(f"u_tau = {u_tau}")
+
+# Generate heights (y)
+heights = np.linspace(height_start, height_end, num_points)
+
+# Calculate u values using smooth wall log law of the wall formula
+u = u_tau * (1. / kappa * np.log(heights * u_tau / kinematic_viscosity) + B)
+print(f"u = {u}")
+
+# Create v values (set all v values to 0)
+v = np.zeros_like(heights)
+w = np.zeros_like(heights)
+mixing_ratio = np.zeros_like(heights)
+potential_temp = np.full_like(heights, 300.0)
+
+# Prepare data to write to file
+data_sounding = np.column_stack((heights, potential_temp, mixing_ratio, u, v))
+data_sponge = np.column_stack((heights, u, v))
+data_inflow = np.column_stack((heights, u, v, w))
+
+# Add the initial rows
+data_sounding = np.vstack([[0.0, 300.0, 0.0, 0.0, 0.0], data_sounding])
+data_sponge = np.vstack([[0.0, 0.0, 0.0], data_sponge])
+data_inflow = np.vstack([[0.0, 0.0, 0.0, 0.0], data_inflow])
+
+# Save data to file
+filename = f"input_ReTau{Re_tau}Ana"
+filetype = ".txt"
+
+np.savetxt(filename + "_sounding" + filetype, data_sounding, fmt="%.8e", delimiter=' ')
+print(f"Data saved to {filename}_sounding{filetype}.")
+np.savetxt(filename + "_sponge" + filetype, data_sponge, fmt="%.8e", delimiter=' ')
+print(f"Data saved to {filename}_sponge{filetype}.")
+np.savetxt(filename + "_inflow" + filetype, data_inflow, fmt="%.8e", delimiter=' ')
+print(f"Data saved to {filename}_inflow{filetype}.")

Exec/DevTests/ABL_perturbation_inflow/incompressible_inputs

@@ -1,37 +1,48 @@
 # ------------------  INPUTS TO MAIN PROGRAM  -------------------
-stop_time = 35.0
-#stop_time = 10.0
+stop_time = 60.0
+#max_step = 5
 
 erf.incompressible = 1
 erf.no_substepping = 1
-
-max_step = 1
-
 amrex.fpe_trap_invalid = 1
-
 fabarray.mfiter_tile_size = 1024 1024 1024
 
 # PROBLEM SIZE & GEOMETRY
-geometry.prob_extent =   2.    0.25    1.
-amr.n_cell           = 128     16     64
+#Larger problem
+geometry.prob_extent =  40     5     10
+amr.n_cell           = 256    32     64
+
+# Quick debug problem
+#geometry.prob_extent =  20     5     10
+#amr.n_cell           =  64    16     32
 
 geometry.is_periodic = 0 1 0
 
-xlo.type = "Outflow"
+xlo.type = "Inflow"
 xhi.type = "Outflow"
-zlo.type = "NoSlipWall"
 zhi.type = "SlipWall"
+#zlo.type = "NoSlipWall"
+
+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
+
+xlo.dirichlet_file = "input_ReTau395Ana_inflow.txt"
+xlo.density =   1.0
+xlo.theta   = 300.0
+#xlo.velocity = 1.0 0.0 0.0
 
 # TIME STEP CONTROL
 erf.cfl = 0.5
-
-erf.dynamicViscosity = 0.0001 # H = 1, therefore utau = 0.2
+erf.dynamicViscosity = 0.001
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
-erf.pert_interval  = 1       # timesteps between perturbation output message
-erf.v              = 1       # verbosity in ERF.cpp
-erf.mg_v           = 1       # verbosity in ERF.cpp
+erf.pert_interval  = 1       # timesteps between perturbation output message XXX
+erf.v              = 0       # verbosity in ERF.cpp XXX
 amr.v              = 0       # verbosity in Amr.cpp
 
 # REFINEMENT / REGRIDDING
@@ -40,48 +51,38 @@ 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      = 10
+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
-
-# Restart for data collection
-#erf.restart = "/home/tma/Desktop/TurbChannel/SpinUp/chk233334"
-#erf.data_log    = "surf.txt" "mean.txt" "flux.txt" "subgrid.txt"
-#erf.profile_int = 33334
+#erf.check_file      = chk    # root name of checkpoint file
+#erf.check_per       = 1.0
 
 # SOLVER CHOICE
 erf.alpha_T = 0.0
-erf.alpha_C = 1.0
-erf.use_gravity = false
-
+erf.alpha_C = 0.0
 erf.molec_diff_type = "None"
 erf.les_type        = "Smagorinsky"
 erf.Cs              = 0.1
-
-erf.buoyancy_type = 1
-
-# Initial condition for the entire field
-erf.init_type = "input_sounding"
-erf.input_sounding_file = "input_ReTau2000DNS_sounding.txt"
+erf.use_gravity     = true
+erf.buoyancy_type   = 2
 
 # Turbulent inflow generation
 erf.perturbation_type = "source"
 erf.perturbation_direction = 1 0 0
 erf.perturbation_layers = 3
-erf.perturbation_offset = 10
+erf.perturbation_offset = 3
 
-erf.perturbation_box_dims = 4 4 8 # 0.25 0.25 0.125
+erf.perturbation_box_dims = 8 8 4
 erf.perturbation_nondimensional = 0.042 # Ri
 erf.perturbation_T_infinity = 300.0
 erf.perturbation_T_intensity = 0.1
 
+# Initial condition for the entire field
+#erf.init_type = "uniform"
+erf.init_type = "input_sounding"
+erf.input_sounding_file = "input_ReTau395Ana_sounding.txt"
+
 # PROBLEM PARAMETERS
 prob.rho_0 = 1.0
-prob.A_0 = 1.0
-prob.U_0 = 0.0
-prob.V_0 = 0.0
-prob.W_0 = 0.0
 prob.T_0 = 300.0