diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst index df7aee683..ad77e8e94 100644 --- 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 index dcca43335..18a91dd96 100644 --- 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 new file mode 100644 index 000000000..2799c7a3d --- /dev/null +++ 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/Exec/DevTests/ABL_perturbation_inflow/toc_inout_inputs b/Exec/DevTests/ABL_perturbation_inflow/compressible_source_inputs similarity index 77% rename from Exec/DevTests/ABL_perturbation_inflow/toc_inout_inputs rename to Exec/DevTests/ABL_perturbation_inflow/compressible_source_inputs index a029e2f28..d6ef78b9e 100644 --- a/Exec/DevTests/ABL_perturbation_inflow/toc_inout_inputs +++ b/Exec/DevTests/ABL_perturbation_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,13 +34,15 @@ 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 @@ -49,7 +50,7 @@ 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" diff --git a/Exec/DevTests/ABL_perturbation_inflow/incompressible_inputs b/Exec/DevTests/ABL_perturbation_inflow/incompressible_direct_inputs similarity index 78% rename from Exec/DevTests/ABL_perturbation_inflow/incompressible_inputs rename to Exec/DevTests/ABL_perturbation_inflow/incompressible_direct_inputs index e7a7abd34..d144bcfc6 100644 --- a/Exec/DevTests/ABL_perturbation_inflow/incompressible_inputs +++ b/Exec/DevTests/ABL_perturbation_inflow/incompressible_direct_inputs @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -stop_time = 60.0 +stop_time = 180.0 #max_step = 5 erf.incompressible = 1 @@ -9,8 +9,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 = 80 2.5 10 +amr.n_cell = 512 16 64 # Quick debug problem #geometry.prob_extent = 20 5 10 @@ -40,8 +40,8 @@ erf.cfl = 0.5 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.sum_interval = 0 # 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 @@ -50,8 +50,8 @@ 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_per_1 = 1.0 +erf.plot_per_1 = 0.1 erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta # CHECKPOINT FILES @@ -68,15 +68,17 @@ erf.use_gravity = true erf.buoyancy_type = 2 # Turbulent inflow generation -erf.perturbation_type = "source" -erf.perturbation_direction = 1 0 0 +#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_offset = 5 -erf.perturbation_box_dims = 8 8 4 +erf.perturbation_box_dims = 16 16 8 erf.perturbation_nondimensional = 0.042 # Ri +#erf.perturbation_nondimensional = 0.1 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" diff --git a/Exec/DevTests/ABL_perturbation_inflow/incompressible_source_inputs b/Exec/DevTests/ABL_perturbation_inflow/incompressible_source_inputs new file mode 100644 index 000000000..c4bd20e4a --- /dev/null +++ b/Exec/DevTests/ABL_perturbation_inflow/incompressible_source_inputs @@ -0,0 +1,90 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +stop_time = 180.0 +#max_step = 5 + +erf.incompressible = 1 +erf.no_substepping = 1 +amrex.fpe_trap_invalid = 1 +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +#Larger problem +geometry.prob_extent = 80 2.5 10 +amr.n_cell = 512 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.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.001 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 0 # timesteps between computing mass +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 +amr.max_level = 0 # maximum level number allowed + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +#erf.plot_per_1 = 1.0 +erf.plot_per_1 = 0.1 +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 = 2 + +# 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 = 5 + +erf.perturbation_box_dims = 8 8 4 +#erf.perturbation_nondimensional = 0.042 # Ri +erf.perturbation_nondimensional = 0.1 +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.T_0 = 300.0 diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 4f9f68a63..a12eeea65 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -249,7 +249,7 @@ struct SolverChoice { amrex::Print() << "Perturb temperature via direct field add\n"; } else if (!turb_pert_type_string.compare("None")) { pert_type = PerturbationType::None; - amrex::Print() << "No perturbation method selected\n"; + amrex::Print() << "No inflow perturbation method selected\n"; } else { amrex::Abort("Dont know this perturbation_type"); } diff --git a/Source/DataStructs/TurbPertStruct.H b/Source/DataStructs/TurbPertStruct.H index 93345f562..2c6884da0 100644 --- a/Source/DataStructs/TurbPertStruct.H +++ b/Source/DataStructs/TurbPertStruct.H @@ -1,10 +1,10 @@ #ifndef _TURB_PERT_STRUCT_H_ #define _TURB_PERT_STRUCT_H_ +#include #include #include #include - /** * Container holding quantities related to turbulent perturbation parameters */ @@ -23,6 +23,9 @@ struct TurbulentPerturbation { // Initializing Perturbation Region // Currently only support perturbation in the x and y direction + // This portion of the code is only called once for initialization, + // therefore efficiency is not required. Dont spend too much time + // here other than to correctly setup the perturbation box region void init_tpi (const int lev, const amrex::IntVect& nx, const amrex::GpuArray dx, @@ -31,6 +34,10 @@ struct TurbulentPerturbation { const int ngrow_state) { + // Initialization for some 0 dependent terms + tpi_Ti = 0.; + tpi_Tinf = 300.; + amrex::ParmParse pp(pp_prefix); // Reading inputs, and placing assertion for the perturbation inflow to work @@ -42,8 +49,6 @@ struct TurbulentPerturbation { pp.query("perturbation_nondimensional",tpi_nonDim); pp.query("perturbation_T_infinity",tpi_Tinf); pp.query("perturbation_T_intensity",tpi_Ti); - pp.query("perturbation_rho_0",tpi_rho); - pp.query("perturbation_cp",tpi_cp); // Check variables message if (tpi_layers < 0) { amrex::Abort("Please provide a valid perturbation layer value (ie. 3-5)"); } @@ -52,42 +57,85 @@ struct TurbulentPerturbation { if (tpi_boxDim[i] == 0) { amrex::Abort("Please provide valid dimensions for perturbation boxes."); } } if (tpi_nonDim < 0.) { amrex::Abort("Please provide a valid nondimensional number (ie. Ri = 0.042)"); } - if (tpi_Tinf < 0.) { amrex::Abort("Please provide a valid ambient temperature value (ie. T_0)"); } - if (tpi_Ti <= 0.) { amrex::Abort("Please provide a valid temperature intensity value (ie. 0-1.0)"); } - if (tpi_rho < 0.) { amrex::Abort("Please provide a valid density (ie. rho_0)"); } - if (tpi_cp < 0.) { amrex::Abort("Please provide a valid specific head capacity"); } - - // Creating perturbation region - amrex::IntVect regionLo; - amrex::IntVect regionHi; - amrex::IntVect boxSize(tpi_boxDim[0],tpi_boxDim[1],tpi_boxDim[2]); // boxSize for individual boxes - - // X-direction perturbation - if (tpi_direction[0]) { - amrex::IntVect tmpLo(0+tpi_offset, 0, 0); - amrex::IntVect tmpHi(((tpi_layers*tpi_boxDim[0])-1)+tpi_offset, nx[1], nx[2]); - regionLo = tmpLo; - regionHi = tmpHi; + if (tpi_Tinf < 0.) { amrex::Abort("Please provide a valid ambient temperature value (ie. T_0 = T_infty)"); } + if (tpi_Ti < 0.) { amrex::Abort("Please provide a valid temperature intensity value (ie. 0-1.0)"); } + + // Creating perturbation regions and initializing with generic size. Temporary size for now + amrex::Box lo_x_bx(amrex::IntVect(0,0,0), amrex::IntVect(1,1,1), amrex::IntVect(0,0,0)); + amrex::Box hi_x_bx(amrex::IntVect(0,0,0), amrex::IntVect(1,1,1), amrex::IntVect(0,0,0)); + amrex::Box lo_y_bx(amrex::IntVect(0,0,0), amrex::IntVect(1,1,1), amrex::IntVect(0,0,0)); + amrex::Box hi_y_bx(amrex::IntVect(0,0,0), amrex::IntVect(1,1,1), amrex::IntVect(0,0,0)); + + // Create a temporary box list to accumulate all the perturbation regions after box modification + amrex::BoxList tmp_bl; + + // boxSize for individual boxes + amrex::IntVect boxSize(tpi_boxDim[0],tpi_boxDim[1],tpi_boxDim[2]); + + // Starting logic to set the size of the perturbation region(s) + amrex::PrintToFile("BoxPerturbationOutput") << "Setting perturbation region in:"; + // ***** X-direction perturbation ***** + if (tpi_direction[0]) { // West + lo_x_bx.setSmall(amrex::IntVect(tpi_offset, tpi_direction[1]*tpi_offset, 0)); + lo_x_bx.setBig(amrex::IntVect((tpi_layers*tpi_boxDim[0]-1)+tpi_offset, nx[1]-(tpi_direction[4]*tpi_offset), nx[2])); + amrex::PrintToFile("BoxPerturbationOutput") << " West face"; + } + + if (tpi_direction[3]) { // East + hi_x_bx.setSmall(amrex::IntVect(nx[0]-((tpi_layers*tpi_boxDim[0]-1)+tpi_offset), tpi_direction[1]*tpi_offset, 0)); + hi_x_bx.setBig(amrex::IntVect(nx[0]-tpi_offset, nx[1]-(tpi_direction[4]*tpi_offset), nx[2])); + amrex::PrintToFile("BoxPerturbationOutput") << " East face"; + } + + // ***** Y-direction Perturbation ***** + if (tpi_direction[1]) { // North + lo_y_bx.setSmall(amrex::IntVect(tpi_direction[0]*tpi_offset, tpi_offset, 0)); + lo_y_bx.setBig(amrex::IntVect(nx[0]-tpi_direction[3]*tpi_offset, ((tpi_layers*tpi_boxDim[1])-1)+tpi_offset, nx[2])); + amrex::PrintToFile("BoxPerturbationOutput") << " North face"; } - // Y-direction Perturbation - if (tpi_direction[1]) { - amrex::IntVect tmpLo(0, 0+tpi_offset, 0); - amrex::IntVect tmpHi(nx[0], ((tpi_layers*tpi_boxDim[1])-1)+tpi_offset, nx[2]); - regionLo = tmpLo; - regionHi = tmpHi; + if (tpi_direction[4]) { // South + hi_y_bx.setSmall(amrex::IntVect(tpi_direction[0]*tpi_offset, nx[1]-((tpi_layers*tpi_boxDim[1]-1)+tpi_offset), 0)); + hi_y_bx.setBig(amrex::IntVect(nx[0]-tpi_direction[3]*tpi_offset, nx[1]-tpi_offset, nx[2])); + amrex::PrintToFile("BoxPerturbationOutput") << " South face"; } - if(tpi_direction[2]) { amrex::Abort("Currently not supporting z-direction flow perturbation"); } + if (tpi_direction[2] || tpi_direction[5]) { amrex::Abort("Currently not supporting z-direction flow perturbation"); } + + // Performing box union for intersecting perturbation regions to avoid overlapping sections (double counting at corners) + if (tpi_direction[0] && tpi_direction[1]) { // Reshaping South smallEnd + amrex::Box lo_x_lo_y_u = lo_x_bx & lo_y_bx; + lo_y_bx.setSmall(amrex::IntVect(lo_x_lo_y_u.bigEnd(0)+1, lo_x_lo_y_u.smallEnd(1), lo_x_lo_y_u.smallEnd(2))); + } + + if (tpi_direction[3] && tpi_direction[1]) { // Reshaping South bigEnd + amrex::Box hi_x_lo_y_u = hi_x_bx & lo_y_bx; + lo_y_bx.setBig(amrex::IntVect(hi_x_lo_y_u.smallEnd(0)-1, hi_x_lo_y_u.bigEnd(1), hi_x_lo_y_u.bigEnd(2))); + } + + if (tpi_direction[0] && tpi_direction[4]) { // Reshaping North smallEnd + amrex::Box lo_x_hi_y_u = lo_x_bx & hi_y_bx; + hi_y_bx.setSmall(amrex::IntVect(lo_x_hi_y_u.bigEnd(0)+1, lo_x_hi_y_u.smallEnd(1), lo_x_hi_y_u.smallEnd(2))); + } + + if (tpi_direction[3] && tpi_direction[4]) { // Reshaping North bigEnd + amrex::Box hi_x_hi_y_u = hi_x_bx & hi_y_bx; + hi_y_bx.setBig(amrex::IntVect(hi_x_hi_y_u.smallEnd(0)-1, hi_x_hi_y_u.bigEnd(1), hi_x_hi_y_u.bigEnd(2))); + } // Creating structure box array for conserved quantity - amrex::Box tmp_bx(regionLo, regionHi, amrex::IntVect(0,0,0)); - amrex::BoxArray tmp_ba(tmp_bx); + if (tpi_direction[0]) { tmp_bl.push_back(lo_x_bx); } + if (tpi_direction[1]) { tmp_bl.push_back(lo_y_bx); } + if (tpi_direction[3]) { tmp_bl.push_back(hi_x_bx); } + if (tpi_direction[4]) { tmp_bl.push_back(hi_y_bx); } + amrex::PrintToFile("BoxPerturbationOutput") << "\nBoxList: " << tmp_bl << "\n"; + amrex::BoxArray tmp_ba(tmp_bl); tmp_ba.maxSize(boxSize); pb_ba.push_back(tmp_ba); // Initializing mean magnitude vector pb_mag.resize(pb_ba[lev].size(), 0.); + pb_netZero.resize(pb_ba[lev].size(), 0.); // Set size of vector and initialize pb_interval.resize(pb_ba[lev].size(), -1.0); @@ -104,6 +152,9 @@ struct TurbulentPerturbation { tpi_Hpb = tpi_boxDim[2]*dx[2]; tpi_lref = sqrt(tpi_Lpb*tpi_Lpb + tpi_Wpb*tpi_Wpb); + tpi_pert_adjust = 0.; + tpi_net_buoyant = 0.; + // Function check point message amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_box_dims: " << tpi_boxDim[0] << " " << tpi_boxDim[1] << " " @@ -111,30 +162,30 @@ struct TurbulentPerturbation { amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_direction: " << tpi_direction[0] << " " << tpi_direction[1] << " " << tpi_direction[2] << "\n\n"; - amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_layers: " << tpi_layers << "\n"; amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_offset: " << tpi_offset << "\n\n"; - amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_nondimensional: " << tpi_nonDim << "\n"; amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_infinity: " << tpi_Tinf << "\n"; amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_T_intensity: " << tpi_Ti << "\n"; - amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_rho_0: " << tpi_rho << "\n"; - amrex::PrintToFile("BoxPerturbationOutput") << "perturbation_cp: " << tpi_cp << "\n\n"; - amrex::PrintToFile("BoxPerturbationOutput") << "Reference length per box = " << tpi_lref << "\n\n"; - amrex::PrintToFile("BoxPerturbationOutput") << "Turbulent perturbation BoxArray:\n" << pb_ba[lev] << "\n"; } + // Perturbation update frequency check. // This function will trigger calc_tpi_meanMag() and calc_tpi_amp(), when - // the local elapsed time is greater than the perturbation frequency + // the local elapsed time is greater than the perturbation frequency. + // At each timestep a total buoyant force sanity check is made to ensure + // that the sum of the buoyant force introduced into the system is net-zero void calc_tpi_update (const int lev, const amrex::Real dt, amrex::MultiFab& mf_xvel, amrex::MultiFab& mf_yvel, amrex::MultiFab& mf_cons) { + // Resettubg the net buoyant force value + tpi_net_buoyant = 0.; + // Setting random number generator for update interval srand( (unsigned) time(NULL) ); @@ -145,12 +196,11 @@ struct TurbulentPerturbation { // Check if the local elapsed time is greater than the update interval if ( pb_interval[boxIdx] <= pb_local_etime[boxIdx] ) { - // Compute mean velocity magnitude within perturbation box + // Compute mean velocity of each perturbation box calc_tpi_meanMag_perBox(boxIdx, lev, mf_xvel, mf_yvel); - //calc_tpi_meanMag_firstCell(boxIdx, lev, mf_xvel, mf_yvel); - // Done for parallelism to avoid Inf being stored in array // Only the rank owning the box will be able to access the storage location + // Done for parallelism to avoid Inf being stored in array if (pb_mag[boxIdx] !=0.) { amrex::Real interval = tpi_lref / pb_mag[boxIdx]; pb_interval[boxIdx] = RandomReal(0.9*interval,1.1*interval); // 10% variation @@ -160,22 +210,7 @@ struct TurbulentPerturbation { calc_tpi_amp(boxIdx, pb_interval[boxIdx]); // Trigger random amplitude storage per cell within perturbation box - // Random perturbation is applied as a white noise on the temperature source term - // it becomes colored through the filtering of temperature transport equation from - // the eddy diffusivity. - for (amrex::MFIter mfi(pb_cell,TileNoZ()); mfi.isValid(); ++mfi) { - amrex::Box vbx = mfi.validbox(); - amrex::Box pbx = convert(pb_ba[lev][boxIdx], m_ixtype); - amrex::Box ubx = pbx & vbx; - if (ubx.ok()) { - amrex::Real amp_copy = pb_amp[boxIdx]; - const amrex::Array4 &pert_cell = pb_cell.array(mfi); - ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { - amrex::Real rand_double = amrex::Random(engine); - pert_cell(i,j,k,0) = (rand_double*2.0 - 1.0) * amp_copy; - }); - } - } + pseudoRandomPert(boxIdx, lev, m_ixtype); // Reset local elapsed time pb_local_etime[boxIdx] = 0.; @@ -183,36 +218,21 @@ struct TurbulentPerturbation { // Increase by timestep of level 0 pb_local_etime[boxIdx] += dt; } // if - } // for - } - -//#define ECKERT_FORMULATION -#define BULK_RICHARDSON_FORMULATION - - // Perturbation amplitude calculation - void calc_tpi_amp (const int boxIdx, const amrex::Real interval) - { - amrex::Real Um = pb_mag[boxIdx]; - pb_amp[boxIdx] = 0.; // Safety step - - #ifdef ECKERT_FORMULATION - // Muñoz-Esparza et al. (2015) Eq. 9 - pb_amp[boxIdx] = Um * Um / (tpi_rho * tpi_cp * tpi_nonDim); - #endif - - #ifdef BULK_RICHARDSON_FORMULATION - amrex::Real beta = 1./tpi_Tinf; // Thermal expansion coefficient - // Pseudo Random temperature the ignores scale when mechanically tripping turbulence - amrex::Real g = (tpi_nonDim * Um * Um) / (tpi_Ti * tpi_Hpb); - //amrex::Real g = CONST_GRAV; + // Per iteration operation of net-zero buoyant force check + if (pb_mag[boxIdx] !=0.) netZeroBuoyantAdd(boxIdx, lev); + tpi_net_buoyant += pb_netZero[boxIdx]; + } // for - // Ma and Senocak (2023) Eq. 8, solving for delta phi - if (g != 0) pb_amp[boxIdx] = (tpi_nonDim * Um * Um) / (g * beta * tpi_Hpb); - #endif + // Normalizing the adjustment based on how many boxes there are + // the values within the array is already normalized by the number + // of cells within each box + tpi_pert_adjust = tpi_net_buoyant / (amrex::Real) pb_ba[lev].size(); - // Ma & Senocak (2023) Eq. 7 - pb_amp[boxIdx] /= interval; + // Per iteration operation of net-zero buoyant force adjustment + for (int boxIdx = 0; boxIdx < pb_ba[lev].size(); boxIdx++) { + if (pb_mag[boxIdx] !=0.) netZeroBuoyantAdjust(boxIdx, lev); + } } //#define INDEX_PERTURB @@ -243,7 +263,103 @@ struct TurbulentPerturbation { } } -// TODO: Test the difference between these two + // Perturbation amplitude calculation + void calc_tpi_amp (const int boxIdx, const amrex::Real interval) + { + amrex::Real Um = pb_mag[boxIdx]; + pb_amp[boxIdx] = 0.; // Safety step + + amrex::Real beta = 1./tpi_Tinf; // Thermal expansion coefficient + + // Pseudo Random temperature the ignores scale when mechanically tripping turbulence + amrex::Real g = CONST_GRAV; + if (tpi_Ti > 0.) g = (tpi_nonDim * Um * Um) / (tpi_Ti * tpi_Hpb); + + // Ma and Senocak (2023) Eq. 8, solving for delta phi + pb_amp[boxIdx] = (tpi_nonDim * Um * Um) / (g * beta * tpi_Hpb); + + if (pt_type == 0) { + // Performing this step converts the perturbation proportionality into + // the forcing term + // Ma & Senocak (2023) Eq. 7 + pb_amp[boxIdx] /= interval; + } + } + + // Assigns pseudo-random (ie. white noise) perturbation to a storage cell, this + // value is then held constant for the duration of the update interval and assigned onto + // the source term + void pseudoRandomPert (const int boxIdx, const int lev, const amrex::IndexType m_ixtype) + { + for (amrex::MFIter mfi(pb_cell,TileNoZ()); mfi.isValid(); ++mfi) { + amrex::Box vbx = mfi.validbox(); + amrex::Box pbx = convert(pb_ba[lev][boxIdx], m_ixtype); + amrex::Box ubx = pbx & vbx; + if (ubx.ok()) { + amrex::Real amp_copy = pb_amp[boxIdx]; + const amrex::Array4 &pert_cell = pb_cell.array(mfi); + ParallelForRNG(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + amrex::Real rand_double = amrex::Random(engine); + pert_cell(i,j,k,0) = (rand_double*2.0 - 1.0) * amp_copy; + }); + } + } + } + + // Checks for net-zero buoyant force introduction into the system + void netZeroBuoyantAdd (const int boxIdx, const int lev) + { + // Creating local copy of PB box array and magnitude + const amrex::BoxArray m_pb_ba = pb_ba[lev]; + amrex::Real* m_pb_netZero = get_pb_netZero(); + + // Create device array for summation + amrex::Vector tmp_h(1,0.); + amrex::Gpu::DeviceVector tmp_d(1,0.); + amrex::Real* tmp = tmp_d.data(); + + // Iterates through the cells of each box and sum the white noise perturbation + for (amrex::MFIter mfi(pb_cell, TileNoZ()) ; mfi.isValid(); ++mfi) { + const amrex::Box &vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + if (ubx.ok()) { + const amrex::Array4 &pert_cell = pb_cell.array(mfi); + int npts = ubx.numPts(); + amrex::Real norm = 1.0 / (amrex::Real) npts; + ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] + AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { + amrex::Gpu::deviceReduceSum(&tmp[0], pert_cell(i,j,k)*norm, handler); + }); + amrex::Gpu::copy(amrex::Gpu::deviceToHost, tmp_d.begin(), tmp_d.end(), tmp_h.begin()); + + // Assigning onto storage array + m_pb_netZero[boxIdx] = tmp_h[0]; + } + } + } + + // If it's not a net-zero buoyant force, then adjust all cell by a normalized value + // to achieve this logic + void netZeroBuoyantAdjust (const int boxIdx, const int lev) + { + // Creating local copy of PB box array and magnitude + const amrex::BoxArray m_pb_ba = pb_ba[lev]; + for (amrex::MFIter mfi(pb_cell, TileNoZ()) ; mfi.isValid(); ++mfi) { + const amrex::Box &vbx = mfi.validbox(); + amrex::Box pbx = convert(m_pb_ba[boxIdx], vbx.ixType()); + amrex::Box ubx = pbx & vbx; + if (ubx.ok()) { + const amrex::Real adjust = tpi_pert_adjust; + const amrex::Array4 &pert_cell = pb_cell.array(mfi); + ParallelFor(ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + pert_cell(i,j,k,0) -= adjust; + }); + } + } + } + +// TODO: Test the difference between these two for Source term perturbation //#define USE_SLAB_AVERAGE #define USE_VOLUME_AVERAGE @@ -258,10 +374,12 @@ struct TurbulentPerturbation { // Creating local copy of PB box array and magnitude const amrex::BoxArray m_pb_ba = pb_ba[lev]; amrex::Real* m_pb_mag = get_pb_mag(); + m_pb_mag[boxIdx] = 0.; // Safety step // Storage of averages per PB // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo) // 2=u (slab_hi), 3=v (slab_hi) + amrex::Vector tmp_h(4,0.); amrex::Gpu::DeviceVector tmp_d(4,0.); amrex::Real* tmp = tmp_d.data(); @@ -277,11 +395,11 @@ struct TurbulentPerturbation { #ifdef USE_VOLUME_AVERAGE int npts = ubx.numPts(); + amrex::Real norm = 1.0 / (amrex::Real) npts; ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k), handler); + amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k)*norm, handler); }); - tmp[0] /= (amrex::Real) npts; #endif // USE_VOLUME_AVERAGE #ifdef USE_SLAB_AVERAGE @@ -289,20 +407,20 @@ struct TurbulentPerturbation { amrex::Box ubxSlab_hi = makeSlab(ubx,2,ubx.bigEnd(2)); int npts_lo = ubxSlab_lo.numPts(); int npts_hi = ubxSlab_hi.numPts(); + amrex::Real norm_lo = 1.0 / (amrex::Real) npts_lo; + amrex::Real norm_hi = 1.0 / (amrex::Real) npts_hi; // Average u in the low slab ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k), handler); + amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k)*norm_lo, handler); }); - tmp[0] /= (amrex::Real) npts_lo; // Average u in the high slab ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[2], xvel_arry(i,j,k), handler); + amrex::Gpu::deviceReduceSum(&tmp[2], xvel_arry(i,j,k)*norm_hi, handler); }); - tmp[2] /= (amrex::Real) npts_hi; #endif // USE_SLAB_AVERAGE } // if } // MFIter @@ -319,11 +437,11 @@ struct TurbulentPerturbation { #ifdef USE_VOLUME_AVERAGE int npts = ubx.numPts(); + amrex::Real norm = 1.0 / (amrex::Real) npts; ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k), handler); + amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k)*norm, handler); }); - tmp[1] /= (amrex::Real) npts; #endif // USE_VOLUME_AVERAGE #ifdef USE_SLAB_AVERAGE @@ -331,105 +449,48 @@ struct TurbulentPerturbation { amrex::Box ubxSlab_hi = makeSlab(ubx,2,ubx.bigEnd(2)); int npts_lo = ubxSlab_lo.numPts(); int npts_hi = ubxSlab_hi.numPts(); + amrex::Real norm_lo = 1.0 / (amrex::Real) npts_lo; + amrex::Real norm_hi = 1.0 / (amrex::Real) npts_hi; // Average v in the low slab ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_lo, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k), handler); + amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k)*norm_lo, handler); }); - tmp[1] /= (amrex::Real) npts_lo; // Average v in the high slab ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubxSlab_hi, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[3], yvel_arry(i,j,k), handler); + amrex::Gpu::deviceReduceSum(&tmp[3], yvel_arry(i,j,k)*norm_hi, handler); }); - tmp[3] /= (amrex::Real) npts_hi; #endif // USE_SLAB_AVERAGE } // if } // MFIter - // Computing the average magnitude within PB - #ifdef USE_SLAB_AVERAGE - m_pb_mag[boxIdx] = 0.5*(sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1]) + sqrt(tmp[2]*tmp[2] + tmp[3]*tmp[3])); - #endif + // Copy from device back to host + amrex::Gpu::copy(amrex::Gpu::deviceToHost, tmp_d.begin(), tmp_d.end(), tmp_h.begin()); + // Computing the average magnitude within PB #ifdef USE_VOLUME_AVERAGE - m_pb_mag[boxIdx] = sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1]); + m_pb_mag[boxIdx] = sqrt(tmp_h[0]*tmp_h[0] + tmp_h[1]*tmp_h[1]); #endif - } - // This emulates calc_tpi_meanMag_perBox() but is heavily simplified for - // Eckert number formulation - void calc_tpi_meanMag_firstCell (const int boxIdx, - const int lev, - amrex::MultiFab& mf_xvel, - amrex::MultiFab& mf_yvel) - { - // Creating local copy of PB box array and magnitude - const amrex::BoxArray m_pb_ba = pb_ba[lev]; - amrex::Real* m_pb_mag = get_pb_mag(); - - // Storage of averages per PB - // Index: 0=u (vol/slab_lo), 1=v (vol/slab_lo) - amrex::Gpu::DeviceVector tmp_d(2,0.); - amrex::Real* tmp = tmp_d.data(); - - // Averaging u components - for (amrex::MFIter mfi(mf_xvel, TileNoZ()) ; mfi.isValid(); ++mfi) { - const amrex::Box &vbx = mfi.validbox(); - // This is hardcoded so that it's always on the 1st cell level (May need to generalize for EB) - amrex::Box pbx(amrex::IntVect(m_pb_ba[boxIdx].smallEnd(0),m_pb_ba[boxIdx].smallEnd(1),0), - amrex::IntVect(m_pb_ba[boxIdx].bigEnd(0) ,m_pb_ba[boxIdx].bigEnd(1) ,1), - vbx.ixType()); - amrex::Box ubx = pbx & vbx; - - // Operation over box union - if (ubx.ok()) { - const amrex::Array4 &xvel_arry = mf_xvel.array(mfi); - int npts = ubx.numPts(); - ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] - AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[0], xvel_arry(i,j,k), handler); - }); - tmp[0] /= (amrex::Real) npts; - } // if - } // MFIter - - // Averaging v components - for (amrex::MFIter mfi(mf_yvel, TileNoZ()) ; mfi.isValid(); ++mfi) { - const amrex::Box &vbx = mfi.validbox(); - // This is hardcoded so that it's always on the 1st cell level (May need to generalize for EB) - amrex::Box pbx(amrex::IntVect(m_pb_ba[boxIdx].smallEnd(0),m_pb_ba[boxIdx].smallEnd(1),0), - amrex::IntVect(m_pb_ba[boxIdx].bigEnd(0) ,m_pb_ba[boxIdx].bigEnd(1) ,1), - vbx.ixType()); - amrex::Box ubx = pbx & vbx; - - // Operation over box union - if (ubx.ok()) { - const amrex::Array4 &yvel_arry = mf_yvel.array(mfi); - int npts = ubx.numPts(); - ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), ubx, [=] - AMREX_GPU_DEVICE(int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept { - amrex::Gpu::deviceReduceSum(&tmp[1], yvel_arry(i,j,k), handler); - }); - tmp[1] /= (amrex::Real) npts; - } // if - } // MFIter - - // Computing the average magnitude within 1st z level slab - m_pb_mag[boxIdx] = sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1]); + #ifdef USE_SLAB_AVERAGE + m_pb_mag[boxIdx] = 0.5*(sqrt(tmp_h[0]*tmp_h[0] + tmp_h[1]*tmp_h[1]) + sqrt(tmp_h[2]*tmp_h[2] + tmp_h[3]*tmp_h[3])); + #endif } // Quality of life function definitions amrex::Real* get_pb_mag() { return pb_mag.data(); } + amrex::Real* get_pb_netZero() { return pb_netZero.data(); } // Output debug message into a file void debug (amrex::Real time) { amrex::PrintToFile("BoxPerturbationOutput") << "#################### PB output at time = " << time << " ####################\n"; - + amrex::PrintToFile("BoxPerturbationOutput") << " Using type: " << pt_type << "\n"; + amrex::PrintToFile("BoxPerturbationOutput") << " Net: " << tpi_net_buoyant << " Adjust : " << tpi_pert_adjust << "\n"; for (int i = 0; i < pb_mag.size(); i++) { amrex::PrintToFile("BoxPerturbationOutput") << "[" << i << "] pb_Umag=" << pb_mag[i] @@ -442,6 +503,8 @@ struct TurbulentPerturbation { std::string pp_prefix {"erf"}; + int pt_type; // TODO: May need to pass in solverChoice to replace this + // Public data members amrex::Vector pb_ba; // PB box array amrex::Vector pb_mag; // BP mean magnitude [m/s] @@ -450,16 +513,8 @@ struct TurbulentPerturbation { // This is after random assignment of equation (10) in Ma and Senocak 2023 amrex::MultiFab pb_cell; - private: - // Random number generation between range (used for interval calculation) - amrex::Real RandomReal (const amrex::Real min, const amrex::Real max) - { - amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX; - return min + r * (max - min); - } - // Private data members int tpi_layers; // Number of layers of perturbation boxes int tpi_offset; // Cells to offset the start of the perturbation region @@ -467,25 +522,32 @@ struct TurbulentPerturbation { amrex::Vector tpi_boxDim; // Dimensions of each perturbation box amrex::Vector tpi_direction; // Direction of the peturbation - amrex::Real tpi_nonDim; // Nondimensional Value - // Richardson Formulation + amrex::Real tpi_nonDim; // Richardson number amrex::Real tpi_Ti; // Temperature intensity value amrex::Real tpi_Tinf; // Reference temperature [K] - // Eckert Formulation - amrex::Real tpi_rho; // Reference density - amrex::Real tpi_cp; // Specific heat capaticty - // Physical dimensions amrex::Real tpi_Hpb; // PB height [m] amrex::Real tpi_Lpb; // PB length [m] amrex::Real tpi_Wpb; // PB width [m] amrex::Real tpi_lref; // PB reference length [m] + amrex::Real tpi_net_buoyant; // Perturbation net-zero calculation storage + amrex::Real tpi_pert_adjust; // Perturbation adjust for net-zero per cell adjustment + // Perturbation data arrays amrex::Vector pb_interval; // PB update time [s] amrex::Vector pb_local_etime; // PB local elapsed time [s] - amrex::Vector pb_amp; // PB perturbation amplitude Ri:[K] Ec:[K/s] + amrex::Vector pb_amp; // PB perturbation amplitude Ri:[K] + amrex::Vector pb_netZero; // PB array used for net zero sum calculation + + // Random number generation between range (used for interval calculation) + amrex::Real RandomReal (const amrex::Real min, const amrex::Real max) + { + amrex::Real r = (amrex::Real) rand() / (amrex::Real) RAND_MAX; + return min + r * (max - min); + } + }; #endif diff --git a/Source/Initialization/ERF_init_TurbPert.cpp b/Source/Initialization/ERF_init_TurbPert.cpp index 458e0abd9..129969b70 100644 --- a/Source/Initialization/ERF_init_TurbPert.cpp +++ b/Source/Initialization/ERF_init_TurbPert.cpp @@ -25,10 +25,19 @@ ERF::turbPert_update (const int lev, const Real local_dt) MultiFab::Copy (xvel_data, lev_new[Vars::xvel], 0, 0, 1, lev_new[Vars::xvel].nGrowVect()); MultiFab::Copy (yvel_data, lev_new[Vars::yvel], 0, 0, 1, lev_new[Vars::yvel].nGrowVect()); + // This logic is done once then stored within TurbPertStruct.H + turbPert.pt_type = -1; + if (solverChoice.pert_type == PerturbationType::perturbSource) { + turbPert.pt_type = 0; + } else if (solverChoice.pert_type == PerturbationType::perturbDirect) { + turbPert.pt_type = 1; + } + AMREX_ALWAYS_ASSERT(turbPert.pt_type >= 0); + // Computing perturbation update time turbPert.calc_tpi_update(lev, local_dt, xvel_data, yvel_data, cons_data); - Print() << "Successfully initialized turbulent perturbation update time and amplitude\n"; + Print() << "Successfully initialized turbulent perturbation update time and amplitude with type: "<< turbPert.pt_type <<"\n"; } // Calculate the perturbation region amplitude. diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index aa897e9c1..d495ea5b8 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -46,18 +46,16 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) // If perturbDirect is selected, directly add the computed perturbation // on the conserved field - // XXX Currently gives "Erroneous arithmetic operation" error - /* if (solverChoice.pert_type == PerturbationType::perturbDirect) { auto m_ixtype = S_old.boxArray().ixType(); // Conserved term for (MFIter mfi(S_old,TileNoZ()); mfi.isValid(); ++mfi) { Box bx = mfi.tilebox(); - const Array4 & cell_data = S_old.array(mfi); - turbPert.apply_tpi(lev, bx, RhoTheta_comp, m_ixtype, cell_data); + const Array4 &cell_data = S_old.array(mfi); + const Array4 &pert_cell = turbPert.pb_cell.array(mfi); + turbPert.apply_tpi(lev, bx, RhoTheta_comp, m_ixtype, cell_data, pert_cell); } } - */ } // configure ABLMost params if used MostWall boundary condition