Merge branch 'development' into ActuatorDisk

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