Synthetic turbulence

docs/documentation/case.md

@@ -1173,6 +1173,29 @@ Usage notes:
 
 - These parameters are for NVIDIA Grace-Hopper and similar architectures with hardware-managed unified memory. They allow MFC to run problems larger than GPU memory by paging data between host and device.
 
+### 20. Synthetic Turbulence
+
+| Parameter                          | Type    | Description                                            |
+| ---:                                | :---:   | :---                                                    |
+| `synthetic_turbulence`              | Logical | Enable synthetic turbulence forcing                    |
+| `synth_seed`                        | Integer | Random seed for wave vector generation                 |
+| `synth_n_shells`                    | Integer | Number of energy shells                                |
+| `synth_U_inf`                       | Real    | Advection velocity for the synthetic turbulence field  |
+| `synth_n_waves_per_shell(s)`        | Integer | Number of random wave vectors in shell `s`             |
+| `synth_k_shell(s)`                  | Real    | Wave-number magnitude of shell `s`                     |
+| `synth_amp_shell(s)`                | Real    | Forcing amplitude of shell `s`                         |
+| `num_turbulent_sources`             | Integer | Number of Gaussian forcing zones                       |
+| `turb_pos(i,1[,2,3])`                | Real    | Center of forcing zone `i` (x[,y,z])                   |
+| `synth_L(i,1[,2,3])`                 | Real    | Full extent of forcing zone `i` (x[,y,z])              |
+
+- `synthetic_turbulence` superimposes a divergence-free, time-advected sum of Fourier modes onto the momentum and energy equations to mimic isotropic turbulence entering the domain.
+
+- Each of the `synth_n_shells` energy shells is defined by a wave-number magnitude `synth_k_shell(s)`, a forcing amplitude `synth_amp_shell(s)`, and `synth_n_waves_per_shell(s)` randomly oriented wave vectors drawn using `synth_seed`.
+
+- `synth_U_inf` advects the turbulent field in the x-direction at a constant velocity.
+
+- Forcing is applied only within Gaussian-windowed zones. Each of the `num_turbulent_sources` zones must specify its center `turb_pos(i,:)` and full extents `synth_L(i,:)` for every active dimension (x, and y, z if `n > 0`/`p > 0`).
+
 ## Enumerations
 
 ### Boundary conditions {#boundary-conditions}

examples/2D_synthetic_turbulence/case.py

@@ -0,0 +1,152 @@
+import json
+import math
+
+Mu = 1.84e-05
+gam_a = 1.4
+gam_b = 1.1
+
+free_stream_vel = 0.5
+
+n_shells = 3
+k_shells = [
+    2 * math.pi / 0.6,  # shell 1 — large eddies
+    2 * math.pi / 0.30,  # shell 2 — medium eddies
+    2 * math.pi / 0.15,
+]  # shell 3 — small eddies
+amp_shells = [
+    0.1 * free_stream_vel,  # m/s  ~5 % of U_inf
+    0.07 * free_stream_vel,  # m/s  ~3 %
+    0.05 * free_stream_vel,
+]  # m/s  ~1 %
+waves_per_shell = [4, 8, 12]  # random directions per shell
+
+# Gaussian source zone
+src_x, src_y = 1.0, 3.0  # center
+src_Lx, src_Ly = 1.0, 4.0  # full extents fed to synth_L
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            # axial direction
+            "x_domain%beg": 0.0e00,
+            "x_domain%end": 12.0,
+            # r direction
+            "y_domain%beg": 0.0e00,
+            "y_domain%end": 6.0,
+            "cyl_coord": "F",
+            "m": 500,
+            "n": 250,
+            "p": 0,
+            "dt": 2.0e-3,
+            "t_step_start": 0,
+            "t_step_stop": 20000,  # 3000
+            "t_step_save": 80,  # 10
+            # Simulation Algorithm Parameters
+            # Only one patches are necessary, the air tube
+            "num_patches": 1,
+            # Use the 5 equation model
+            "model_eqns": 2,
+            # 6 equations model does not need the K \div(u) term
+            "alt_soundspeed": "F",
+            # One fluids: air
+            "num_fluids": 2,
+            # time step
+            "mpp_lim": "F",
+            # Correct errors when computing speed of sound
+            "mixture_err": "T",
+            # Use TVD RK3 for time marching
+            "time_stepper": 3,
+            # Reconstruct the primitive variables to minimize spurious
+            # Use WENO5
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "weno_Re_flux": "F",
+            "weno_avg": "T",
+            "avg_state": 2,
+            # Use the mapped WENO weights to maintain monotinicity
+            "mapped_weno": "T",
+            "null_weights": "F",
+            "mp_weno": "T",
+            # Use the HLLC  Riemann solver
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            # We use reflective boundary conditions at octant edges and
+            # non-reflective boundary conditions at the domain edges
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "bc_y%beg": -3,
+            "bc_y%end": -3,
+            # Set IB to True and add 1 patch
+            "ib": "T",
+            "num_ibs": 1,
+            "fd_order": 2,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "E_wrt": "T",
+            "parallel_io": "T",
+            # Patch: Constant Tube filled with air
+            # Specify the cylindrical air tube grid geometry
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 6.0,
+            # Uniform medium density, centroid is at the center of the domain
+            "patch_icpp(1)%y_centroid": 3.0,
+            "patch_icpp(1)%length_x": 12.0,
+            "patch_icpp(1)%length_y": 6.0,
+            # Specify the patch primitive variables
+            "patch_icpp(1)%vel(1)": free_stream_vel,
+            "patch_icpp(1)%vel(2)": 0.0e00,
+            "patch_icpp(1)%pres": 1.0e00,
+            "patch_icpp(1)%alpha_rho(1)": 0.8e00,
+            "patch_icpp(1)%alpha(1)": 0.8e00,
+            "patch_icpp(1)%alpha_rho(2)": 0.2e00,
+            "patch_icpp(1)%alpha(2)": 0.2e00,
+            # Patch: Airfoil Immersed Boundary
+            "patch_ib(1)%geometry": 4,
+            "patch_ib(1)%x_centroid": 6.0,
+            "patch_ib(1)%y_centroid": 3.0,
+            "patch_ib(1)%airfoil_id": 1,
+            "ib_airfoil(1)%c": 2.0,
+            "ib_airfoil(1)%t": 0.15,
+            "ib_airfoil(1)%p": 0.4,
+            "ib_airfoil(1)%m": 0.02,
+            "patch_ib(1)%angles(3)": -0.5235987756,  # 30 degrees clockwise rotation, in radians
+            "patch_ib(1)%angular_vel(3)": "15.0 * 0.1 * pi * sin(0.1 * pi * t) * pi / 180.",
+            "patch_ib(1)%moving_ibm": 1,
+            # Fluids Physical Parameters
+            # Use the same stiffness as the air bubble
+            "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50 (Not 1.40)
+            "fluid_pp(1)%pi_inf": 0,
+            "fluid_pp(2)%gamma": 1.0e00 / (gam_b - 1.0e00),  # 2.50 (Not 1.40)
+            "fluid_pp(2)%pi_inf": 0,
+            # -- Synthetic turbulence --
+            "synthetic_turbulence": "T",
+            "synth_seed": 42,
+            "synth_n_shells": n_shells,
+            "synth_U_inf": free_stream_vel,
+            "num_turbulent_sources": 1,
+            # Energy shells (wave-number magnitudes, amplitudes, wave counts)
+            "synth_k_shell(1)": k_shells[0],
+            "synth_k_shell(2)": k_shells[1],
+            "synth_k_shell(3)": k_shells[2],
+            "synth_amp_shell(1)": amp_shells[0],
+            "synth_amp_shell(2)": amp_shells[1],
+            "synth_amp_shell(3)": amp_shells[2],
+            "synth_n_waves_per_shell(1)": waves_per_shell[0],
+            "synth_n_waves_per_shell(2)": waves_per_shell[1],
+            "synth_n_waves_per_shell(3)": waves_per_shell[2],
+            # Gaussian forcing zone (source 1)
+            "turb_pos(1,1)": src_x,
+            "turb_pos(1,2)": src_y,
+            "turb_pos(1,3)": 0.0,  # z not used in 2-D
+            "synth_L(1,1)": src_Lx,
+            "synth_L(1,2)": src_Ly,
+            "synth_L(1,3)": 1.0,  # z extent irrelevant in 2-D; set nonzero to avoid div-by-zero
+        }
+    )
+)

src/common/m_constants.fpp

@@ -114,6 +114,10 @@ module m_constants
     integer, parameter :: BC_NO_SLIP_WALL = -16
     integer, parameter :: BC_DIRICHLET = -17
 
+    ! Synthetic turbulence array size limits
+    integer, parameter :: num_synth_shells_max = 50  !< Max energy shells for synthetic turbulence
+    integer, parameter :: num_turb_sources_max = 10  !< Max Gaussian forcing zones for synthetic turbulence
+
     ! Named values for enumerated case parameters (e.g. riemann_solver_hllc).
     ! AUTO-GENERATED from "names" in toolchain/mfc/params/definitions.py.
     #:include 'generated_constants.fpp'

src/common/m_global_parameters_common.fpp

@@ -98,6 +98,7 @@ module m_global_parameters_common
     $:GPU_DECLARE(create='[Bx0]')
     $:GPU_DECLARE(create='[tau_star, cont_damage_s, alpha_bar]')
     $:GPU_DECLARE(create='[hyper_cleaning_speed, hyper_cleaning_tau]')
+    $:GPU_DECLARE(create='[synthetic_turbulence, num_turbulent_sources, synth_U_inf]')
     #:if not MFC_CASE_OPTIMIZATION
         $:GPU_DECLARE(create='[num_dims, num_vels, weno_polyn, weno_order]')
         $:GPU_DECLARE(create='[weno_num_stencils, num_fluids, wenojs]')

src/common/m_helper.fpp

@@ -18,7 +18,7 @@ module m_helper
     public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, &
         & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, &
         & f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, &
-        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product
+        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul
 
 contains
 
@@ -297,6 +297,48 @@ contains
 
     end function f_cross
 
+    !> Generate a unit vector uniformly distributed on the sphere from two random parameters.
+    function f_unit_vector(theta, eta) result(vec)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in)   :: theta, eta
+        real(wp)               :: zeta, xi
+        real(wp), dimension(3) :: vec
+
+        xi = 2._wp*pi*theta
+        zeta = acos(2._wp*eta - 1._wp)
+        vec(1) = sin(zeta)*cos(xi)
+        vec(2) = sin(zeta)*sin(xi)
+        vec(3) = cos(zeta)
+
+    end function f_unit_vector
+
+    !> Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
+    subroutine s_prng(var, seed)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        integer, intent(inout) :: seed
+        real(wp), intent(out)  :: var
+
+        seed = mod(modmul(seed), modulus)
+        var = seed/real(modulus, wp)
+
+    end subroutine s_prng
+
+    !> Compute a modular multiplication step for the linear congruential pseudo-random number generator.
+    function modmul(a) result(val)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        integer, intent(in) :: a
+        integer             :: val
+        real(wp)            :: x, y
+
+        x = (multiplier/real(modulus, wp))*a + (increment/real(modulus, wp))
+        y = nint((x - floor(x))*decimal_trim)/decimal_trim
+        val = nint(y*modulus)
+
+    end function modmul
+
     !> @brief Computes the cross product c = a x b of two 3D vectors.
     subroutine s_cross_product(a, b, c)
 

src/pre_process/m_perturbation.fpp

@@ -334,45 +334,6 @@ contains
 
     end subroutine s_generate_random_perturbation
 
-    !> Generate a unit vector uniformly distributed on the sphere from two random parameters.
-    function f_unit_vector(theta, eta) result(vec)
-
-        real(wp), intent(in)   :: theta, eta
-        real(wp)               :: zeta, xi
-        real(wp), dimension(3) :: vec
-
-        xi = 2._wp*pi*theta
-        zeta = acos(2._wp*eta - 1._wp)
-        vec(1) = sin(zeta)*cos(xi)
-        vec(2) = sin(zeta)*sin(xi)
-        vec(3) = cos(zeta)
-
-    end function f_unit_vector
-
-    !> Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
-    subroutine s_prng(var, seed)
-
-        integer, intent(inout) :: seed
-        real(wp), intent(out)  :: var
-
-        seed = mod(modmul(seed), modulus)
-        var = seed/real(modulus, wp)
-
-    end subroutine s_prng
-
-    !> Compute a modular multiplication step for the linear congruential pseudo-random number generator.
-    function modmul(a) result(val)
-
-        integer, intent(in) :: a
-        integer             :: val
-        real(wp)            :: x, y
-
-        x = (multiplier/real(modulus, wp))*a + (increment/real(modulus, wp))
-        y = nint((x - floor(x))*decimal_trim)/decimal_trim
-        val = nint(y*modulus)
-
-    end function modmul
-
     !> Deallocate the temporary primitive variable array used by elliptic smoothing.
     impure subroutine s_finalize_perturbation_module()