diff --git a/docs/documentation/case.md b/docs/documentation/case.md index b62b3adcbd..fd162fe8af 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -1111,6 +1111,29 @@ This boundary condition can be used for fixed-temperature (isothermal) walls at - 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} diff --git a/examples/2D_synthetic_turbulence/case.py b/examples/2D_synthetic_turbulence/case.py new file mode 100644 index 0000000000..1e2b3d9f4c --- /dev/null +++ b/examples/2D_synthetic_turbulence/case.py @@ -0,0 +1,151 @@ +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, + # 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 + } + ) +) diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 5452c13a21..8a41fda7c1 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -117,6 +117,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' diff --git a/src/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp index 1f60d2cae3..7b35c3fbda 100644 --- a/src/simulation/m_body_forces.fpp +++ b/src/simulation/m_body_forces.fpp @@ -10,6 +10,7 @@ module m_body_forces use m_derived_types use m_global_parameters use m_variables_conversion + use m_mpi_proxy use m_nvtx ! $:USE_GPU_MODULE() @@ -17,16 +18,37 @@ module m_body_forces implicit none private - public :: s_compute_body_forces_rhs, s_initialize_body_forces_module, s_finalize_body_forces_module + public :: s_compute_body_forces_rhs, s_compute_synthetic_forces_rhs, s_initialize_body_forces_module, & + & s_finalize_body_forces_module real(wp), allocatable, dimension(:,:,:) :: rhoM $:GPU_DECLARE(create='[rhoM]') + integer :: num_synthetic_wave_numbers + $:GPU_DECLARE(create='[num_synthetic_wave_numbers]') + + real(wp), allocatable, dimension(:) :: synthetic_k_x, synthetic_k_y, synthetic_k_z + real(wp), allocatable, dimension(:) :: synthetic_phase, synthetic_amp + ! Solenoidal (divergence-free) unit vector per mode: perpendicular to k in the plane of excitation. + ! 1-D: ex = 1, ey = 0 (streamwise only) + ! 2-D: ex = -k_y/|k|, ey = k_x/|k| (perpendicular to k on the circle) + ! 3-D: stored from random frame built at init + real(wp), allocatable, dimension(:) :: synthetic_ex, synthetic_ey, synthetic_ez + $:GPU_DECLARE(create='[synthetic_k_x, synthetic_k_y, synthetic_k_z, synthetic_phase, synthetic_amp]') + $:GPU_DECLARE(create='[synthetic_ex, synthetic_ey, synthetic_ez]') + contains - !> Initialize the body forces module + !> Initialize the body forces module. When synthetic_turbulence is enabled, + !> generates random wave vectors for each energy shell on rank 0, then + !> broadcasts to all MPI ranks and copies to GPU. impure subroutine s_initialize_body_forces_module + integer :: s, m_wave, m_global + integer :: n_seed + integer, allocatable :: seed_arr(:) + real(wp) :: rn1, rn2, kz, r_xy, k_mag + if (n > 0) then if (p > 0) then @:ALLOCATE(rhoM(-buff_size:buff_size + m, -buff_size:buff_size + n, -buff_size:buff_size + p)) @@ -37,6 +59,122 @@ contains @:ALLOCATE(rhoM(-buff_size:buff_size + m, 0:0, 0:0)) end if + if (.not. synthetic_turbulence) return + if (synth_n_shells <= 0) return + + num_synthetic_wave_numbers = sum(synth_n_waves_per_shell(1:synth_n_shells)) + + if (num_synthetic_wave_numbers == 0) return + + @:ALLOCATE(synthetic_k_x(1:num_synthetic_wave_numbers)) + @:ALLOCATE(synthetic_phase(1:num_synthetic_wave_numbers)) + @:ALLOCATE(synthetic_amp(1:num_synthetic_wave_numbers)) + @:ALLOCATE(synthetic_ex(1:num_synthetic_wave_numbers)) + @:ALLOCATE(synthetic_ey(1:num_synthetic_wave_numbers)) + @:ALLOCATE(synthetic_ez(1:num_synthetic_wave_numbers)) + + if (num_dims > 1) then + @:ALLOCATE(synthetic_k_y(1:num_synthetic_wave_numbers)) + end if + + if (num_dims == 3) then + @:ALLOCATE(synthetic_k_z(1:num_synthetic_wave_numbers)) + end if + + ! Generate random wave vectors and phases on rank 0, then broadcast + if (proc_rank == 0) then + call random_seed(size=n_seed) + allocate (seed_arr(n_seed)) + ! Vary each element so the seed state is not degenerate + do s = 1, n_seed + seed_arr(s) = synth_seed + (s - 1)*1000003 + end do + call random_seed(put=seed_arr) + deallocate (seed_arr) + + m_global = 0 + do s = 1, synth_n_shells + k_mag = synth_k_shell(s) + do m_wave = 1, synth_n_waves_per_shell(s) + m_global = m_global + 1 + + if (num_dims == 1) then + synthetic_k_x(m_global) = k_mag + ! Streamwise forcing only in 1-D + synthetic_ex(m_global) = 1._wp + synthetic_ey(m_global) = 0._wp + synthetic_ez(m_global) = 0._wp + else if (num_dims == 2) then + ! Azimuthal angle theta uniform on [0, 2*pi) + call random_number(rn1) + rn1 = rn1*2._wp*pi + synthetic_k_x(m_global) = k_mag*cos(rn1) + synthetic_k_y(m_global) = k_mag*sin(rn1) + ! Solenoidal direction perpendicular to k: (-sin theta, cos theta) + synthetic_ex(m_global) = -sin(rn1) + synthetic_ey(m_global) = cos(rn1) + synthetic_ez(m_global) = 0._wp + else + ! Uniform on sphere: azimuthal phi, cos(polar) uniform in [-1,1] + call random_number(rn1) + call random_number(rn2) + rn1 = rn1*2._wp*pi + kz = 2._wp*rn2 - 1._wp + r_xy = sqrt(max(0._wp, 1._wp - kz*kz)) + synthetic_k_x(m_global) = k_mag*r_xy*cos(rn1) + synthetic_k_y(m_global) = k_mag*r_xy*sin(rn1) + synthetic_k_z(m_global) = k_mag*kz + ! Solenoidal direction: random unit vector perpendicular to k. + ! Build an orthonormal frame: k_hat x (0,0,1) if not degenerate. + if (abs(kz) < 0.9_wp) then + ! Cross k_hat with z-hat -> (-k_y, k_x, 0) / r_xy + synthetic_ex(m_global) = -sin(rn1) + synthetic_ey(m_global) = cos(rn1) + synthetic_ez(m_global) = 0._wp + else + ! k nearly parallel to z; cross with x-hat instead k_hat x x_hat = (0, -kz, k_y)/norm (normalised) + synthetic_ex(m_global) = 0._wp + synthetic_ey(m_global) = -kz/max(sqrt(kz*kz + (k_mag*sin(rn1))**2), 1e-10_wp) + synthetic_ez(m_global) = (k_mag*sin(rn1))/max(sqrt(kz*kz + (k_mag*sin(rn1))**2), 1e-10_wp) + end if + end if + + call random_number(rn1) + synthetic_phase(m_global) = rn1*2._wp*pi + + synthetic_amp(m_global) = synth_amp_shell(s) + end do + end do + end if + + ! Broadcast from rank 0 to all ranks + call s_mpi_send_random_number(synthetic_k_x, num_synthetic_wave_numbers) + call s_mpi_send_random_number(synthetic_phase, num_synthetic_wave_numbers) + call s_mpi_send_random_number(synthetic_amp, num_synthetic_wave_numbers) + call s_mpi_send_random_number(synthetic_ex, num_synthetic_wave_numbers) + call s_mpi_send_random_number(synthetic_ey, num_synthetic_wave_numbers) + call s_mpi_send_random_number(synthetic_ez, num_synthetic_wave_numbers) + + if (num_dims > 1) then + call s_mpi_send_random_number(synthetic_k_y, num_synthetic_wave_numbers) + end if + + if (num_dims == 3) then + call s_mpi_send_random_number(synthetic_k_z, num_synthetic_wave_numbers) + end if + + ! Push wave data to GPU + $:GPU_UPDATE(device='[num_synthetic_wave_numbers, synthetic_k_x, synthetic_phase, synthetic_amp]') + $:GPU_UPDATE(device='[synthetic_ex, synthetic_ey, synthetic_ez]') + + if (num_dims > 1) then + $:GPU_UPDATE(device='[synthetic_k_y]') + end if + + if (num_dims == 3) then + $:GPU_UPDATE(device='[synthetic_k_z]') + end if + end subroutine s_initialize_body_forces_module !> Compute the acceleration at time t @@ -58,7 +196,7 @@ contains subroutine s_compute_mixture_density(q_cons_vf) type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf - integer :: i, j, k, l !< standard iterators + integer :: i, j, k, l $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3) do l = 0, p @@ -81,7 +219,7 @@ contains type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - integer :: i, j, k, l !< Loop variables + integer :: i, j, k, l call s_compute_acceleration(mytime) call s_compute_mixture_density(q_cons_vf) @@ -98,8 +236,7 @@ contains end do $:END_GPU_PARALLEL_LOOP() - if (bf_x) then ! x-direction body forces - + if (bf_x) then $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3) do l = 0, p do k = 0, n @@ -113,8 +250,7 @@ contains $:END_GPU_PARALLEL_LOOP() end if - if (bf_y) then ! y-direction body forces - + if (bf_y) then $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3) do l = 0, p do k = 0, n @@ -129,8 +265,7 @@ contains $:END_GPU_PARALLEL_LOOP() end if - if (bf_z) then ! z-direction body forces - + if (bf_z) then $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3) do l = 0, p do k = 0, n @@ -146,11 +281,138 @@ contains end subroutine s_compute_body_forces_rhs + !> Compute the synthetic turbulence RHS contribution. + !> + !> Each Fourier mode m has wave vector k_m and a pre-computed solenoidal + !> (divergence-free) direction e_m perpendicular to k_m. The force field + !> is F = rho * (sum_m A_m * cos(k_m.(x-U*t*xhat) + phi_m) * e_m) * (U/Lx) * G. + !> This excites all velocity components correctly and is zero outside the + !> bounding box of the Gaussian window. + subroutine s_compute_synthetic_forces_rhs(q_prim_vf, q_cons_vf, rhs_vf) + + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf + type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf + integer :: i, j, k, l, n_mode, turb_idx + real(wp) :: pos_x, pos_y, pos_z + real(wp) :: gauss_env, G_norm, f_scale + real(wp) :: a_m, phase_arg + real(wp) :: force_x, force_y, force_z + real(wp) :: rho_local, u_local_x, u_local_y, u_local_z + real(wp) :: adv_offset + logical :: in_box + + ! G_bar for exp(-pi/2 * (2x/L)^2) over [-L/2, L/2]: erf(sqrt(pi/2))/sqrt(pi/2) ~ 0.7602 + real(wp), parameter :: G_BAR = 0.760173_wp + + ! Zero momentum and energy RHS for the synthetic turbulence pass + + $:GPU_PARALLEL_LOOP(private='[i, j, k, l]', collapse=4) + do i = eqn_idx%mom%beg, eqn_idx%E + do l = 0, p + do k = 0, n + do j = 0, m + rhs_vf(i)%sf(j, k, l) = 0._wp + end do + end do + end do + end do + $:END_GPU_PARALLEL_LOOP() + + ! Pre-compute advection offset on CPU (mytime has no GPU declaration) + adv_offset = synth_U_inf*mytime + + do turb_idx = 1, num_turbulent_sources + $:GPU_PARALLEL_LOOP(private='[j, k, l, n_mode, pos_x, pos_y, pos_z, gauss_env, G_norm, f_scale, a_m, phase_arg, & + & force_x, force_y, force_z, rho_local, u_local_x, u_local_y, u_local_z, in_box]', collapse=3) + do l = 0, p + do k = 0, n + do j = 0, m + ! Position relative to forcing-zone center + pos_x = x_cc(j) - turb_pos(turb_idx, 1) + pos_y = 0._wp + pos_z = 0._wp + if (num_dims > 1) pos_y = y_cc(k) - turb_pos(turb_idx, 2) + if (num_dims == 3) pos_z = z_cc(l) - turb_pos(turb_idx, 3) + + ! Bounding-box skip + in_box = abs(pos_x) <= 0.5_wp*synth_L(turb_idx, 1) + if (num_dims > 1) in_box = in_box .and. abs(pos_y) <= 0.5_wp*synth_L(turb_idx, 2) + if (num_dims == 3) in_box = in_box .and. abs(pos_z) <= 0.5_wp*synth_L(turb_idx, 3) + + if (in_box) then + ! Gaussian envelope, normalised by G_BAR^num_dims + gauss_env = exp(-0.5_wp*pi*(2._wp*pos_x/synth_L(turb_idx, 1))**2) + if (num_dims > 1) gauss_env = gauss_env*exp(-0.5_wp*pi*(2._wp*pos_y/synth_L(turb_idx, 2))**2) + if (num_dims == 3) gauss_env = gauss_env*exp(-0.5_wp*pi*(2._wp*pos_z/synth_L(turb_idx, 3))**2) + G_norm = gauss_env/G_BAR**num_dims + + ! Accumulate solenoidal force components over all modes. + ! Each mode contributes A_m * cos(phase) in direction e_m + ! (e_m is perpendicular to k_m, precomputed at init). + force_x = 0._wp + force_y = 0._wp + force_z = 0._wp + do n_mode = 1, num_synthetic_wave_numbers + phase_arg = synthetic_k_x(n_mode)*(x_cc(j) - adv_offset) + synthetic_phase(n_mode) + if (num_dims > 1) phase_arg = phase_arg + synthetic_k_y(n_mode)*y_cc(k) + if (num_dims == 3) phase_arg = phase_arg + synthetic_k_z(n_mode)*z_cc(l) + a_m = synthetic_amp(n_mode)*cos(phase_arg) + force_x = force_x + a_m*synthetic_ex(n_mode) + force_y = force_y + a_m*synthetic_ey(n_mode) + force_z = force_z + a_m*synthetic_ez(n_mode) + end do + + ! Scale: rho * (U_inf/L_x) * G_norm + rho_local = q_prim_vf(eqn_idx%cont%beg)%sf(j, k, l) + f_scale = rho_local*(synth_U_inf/synth_L(turb_idx, 1))*G_norm + + force_x = force_x*f_scale + force_y = force_y*f_scale + force_z = force_z*f_scale + + ! Local velocities for energy update F.u + u_local_x = q_prim_vf(eqn_idx%mom%beg)%sf(j, k, l) + u_local_y = 0._wp + u_local_z = 0._wp + if (num_dims > 1) u_local_y = q_prim_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) + if (num_dims == 3) u_local_z = q_prim_vf(eqn_idx%mom%end)%sf(j, k, l) + + rhs_vf(eqn_idx%mom%beg)%sf(j, k, l) = rhs_vf(eqn_idx%mom%beg)%sf(j, k, l) + force_x + if (num_dims > 1) rhs_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = rhs_vf(eqn_idx%mom%beg + 1)%sf(j, k, & + & l) + force_y + if (num_dims == 3) rhs_vf(eqn_idx%mom%end)%sf(j, k, l) = rhs_vf(eqn_idx%mom%end)%sf(j, k, l) + force_z + rhs_vf(eqn_idx%E)%sf(j, k, l) = rhs_vf(eqn_idx%E)%sf(j, k, & + & l) + force_x*u_local_x + force_y*u_local_y + force_z*u_local_z + end if + end do + end do + end do + $:END_GPU_PARALLEL_LOOP() + end do + + end subroutine s_compute_synthetic_forces_rhs + !> Finalize the body forces module impure subroutine s_finalize_body_forces_module @:DEALLOCATE(rhoM) + if (num_synthetic_wave_numbers > 0) then + @:DEALLOCATE(synthetic_k_x) + @:DEALLOCATE(synthetic_phase) + @:DEALLOCATE(synthetic_amp) + @:DEALLOCATE(synthetic_ex) + @:DEALLOCATE(synthetic_ey) + @:DEALLOCATE(synthetic_ez) + if (num_dims > 1) then + @:DEALLOCATE(synthetic_k_y) + end if + if (num_dims == 3) then + @:DEALLOCATE(synthetic_k_z) + end if + end if + end subroutine s_finalize_body_forces_module end module m_body_forces diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index 2ed8aca3c3..f3e733ceb0 100644 --- a/src/simulation/m_checker.fpp +++ b/src/simulation/m_checker.fpp @@ -42,6 +42,10 @@ contains call s_check_inputs_particle_clouds end if + if (synthetic_turbulence) then + call s_check_inputs_synthetic_turbulence + end if + end subroutine s_check_inputs !> Checks constraints on compiler options @@ -127,4 +131,26 @@ contains end subroutine s_check_inputs_particle_clouds + !> Checks that each active synthetic-turbulence forcing zone has a fully specified position and a positive size in every active + !! dimension + impure subroutine s_check_inputs_synthetic_turbulence + + integer :: i, d + character(len=5) :: idxStr + + @:PROHIBIT(num_turbulent_sources <= 0, "num_turbulent_sources must be > 0 when synthetic_turbulence is enabled") + + do i = 1, num_turbulent_sources + call s_int_to_str(i, idxStr) + do d = 1, num_dims + @:PROHIBIT(f_is_default(turb_pos(i, d)), & + & "turb_pos("//trim(idxStr) & + & //",:) must be specified for all num_dims when synthetic_turbulence is enabled") + @:PROHIBIT(f_is_default(synth_L(i, d)) .or. synth_L(i, d) <= 0._wp, & + & "synth_L("//trim(idxStr)//",:) must be positive for all num_dims when synthetic_turbulence is enabled") + end do + end do + + end subroutine s_check_inputs_synthetic_turbulence + end module m_checker diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 77eb90a31d..d623076847 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -128,6 +128,17 @@ module m_global_parameters $:GPU_DECLARE(create='[accel_bf]') ! $:GPU_DECLARE(create='[k_x,w_x,p_x,g_x,k_y,w_y,p_y,g_y,k_z,w_z,p_z,g_z]') + ! Synthetic turbulence + logical :: synthetic_turbulence + integer :: synth_seed, synth_n_shells, num_turbulent_sources + real(wp) :: synth_U_inf + integer, dimension(num_synth_shells_max) :: synth_n_waves_per_shell + real(wp), dimension(num_synth_shells_max) :: synth_k_shell, synth_amp_shell + real(wp), dimension(num_turb_sources_max, 3) :: turb_pos, synth_L + $:GPU_DECLARE(create='[synthetic_turbulence, num_turbulent_sources]') + $:GPU_DECLARE(create='[synth_U_inf, synth_n_waves_per_shell, synth_k_shell, synth_amp_shell]') + $:GPU_DECLARE(create='[turb_pos, synth_L]') + integer :: cpu_start, cpu_end, cpu_rate #:if not MFC_CASE_OPTIMIZATION @@ -351,7 +362,7 @@ module m_global_parameters type(pres_field), allocatable, dimension(:) :: pb_ts type(pres_field), allocatable, dimension(:) :: mv_ts - $:GPU_DECLARE(create='[pb_ts, mv_ts]') + $:GPU_DECLARE(create='[mytime, pb_ts, mv_ts]') !> @name lagrangian subgrid bubble parameters !> @{! @@ -586,6 +597,17 @@ contains #:endfor #:endfor + synthetic_turbulence = .false. + synth_seed = 1234 + synth_n_shells = dflt_int + num_turbulent_sources = 0 + synth_U_inf = dflt_real + synth_n_waves_per_shell = 0 + synth_k_shell = dflt_real + synth_amp_shell = dflt_real + turb_pos = dflt_real + synth_L = dflt_real + fft_wrt = .false. do j = 1, num_probes_max @@ -1152,6 +1174,12 @@ contains $:GPU_UPDATE(device='[relax, relax_model, palpha_eps, ptgalpha_eps]') + if (synthetic_turbulence) then + $:GPU_UPDATE(device='[synthetic_turbulence, num_turbulent_sources]') + $:GPU_UPDATE(device='[synth_U_inf, synth_n_waves_per_shell, synth_k_shell, synth_amp_shell]') + $:GPU_UPDATE(device='[turb_pos, synth_L]') + end if + ! Allocating grid variables for the x-, y- and z-directions @:ALLOCATE(x_cb(-1 - buff_size:m + buff_size)) @:ALLOCATE(x_cc(-buff_size:m + buff_size)) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 99aeb1e382..76c6b872a4 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1383,6 +1383,10 @@ contains local_ib_idx_old(i) = patch_ib(local_ib_patch_ids(i))%gbl_patch_id end do + ! Sync GPU-updated fields (angles, angular_vel, centroids) to host before + ! compaction and MPI packing, which read from host memory. + $:GPU_UPDATE(host='[patch_ib]') + ! delete any particles that no longer need to be tracked and coalesce the array output_idx = 0 local_output_idx = 0 diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index ae955cef55..f83f454bde 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -96,7 +96,7 @@ contains & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'surface_tension', & & 'shear_stress', 'bulk_stress', 'bubbles_lagrange', & & 'hyperelasticity', 'down_sample', 'fft_wrt', & - & 'hyper_cleaning', 'ib_state_wrt'] + & 'hyper_cleaning', 'ib_state_wrt', 'synthetic_turbulence'] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -265,6 +265,17 @@ contains call MPI_BCAST(nv_uvm_out_of_core, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(nv_uvm_igr_temps_on_gpu, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(nv_uvm_pref_gpu, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + + ! Synthetic turbulence + call MPI_BCAST(synth_seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(synth_n_shells, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(num_turbulent_sources, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(synth_U_inf, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(synth_n_waves_per_shell, num_synth_shells_max, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(synth_k_shell, num_synth_shells_max, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(synth_amp_shell, num_synth_shells_max, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(turb_pos, num_turb_sources_max*3, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(synth_L, num_turb_sources_max*3, mpi_p, 0, MPI_COMM_WORLD, ierr) #endif end subroutine s_mpi_bcast_user_inputs diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 7a4a1a7df8..043b3ed687 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -895,7 +895,7 @@ contains call s_write_ib_state_file(0) end if end if - if (bodyForces) call s_initialize_body_forces_module() + if (bodyForces .or. synthetic_turbulence) call s_initialize_body_forces_module() if (acoustic_source) call s_precalculate_acoustic_spatial_sources() ! Initialize the Temperature cache. @@ -1101,7 +1101,7 @@ contains call s_finalize_mpi_proxy_module() if (surface_tension) call s_finalize_surface_tension_module() - if (bodyForces) call s_finalize_body_forces_module() + if (bodyForces .or. synthetic_turbulence) call s_finalize_body_forces_module() if (ib) call s_finalize_ibm_module() call s_mpi_finalize() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 9a59f955f2..89fa4d8c45 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -539,8 +539,12 @@ contains $:END_GPU_PARALLEL_LOOP() end if + $:GPU_UPDATE(device='[mytime]') if (bodyForces) call s_apply_bodyforces(q_cons_ts(1)%vf, q_prim_vf, rhs_vf, rk_coef(s, 3)*dt/rk_coef(s, 4)) + if (synthetic_turbulence) call s_apply_synthetic_turbulence_force(q_cons_ts(1)%vf, q_prim_vf, rhs_vf, rk_coef(s, & + & 3)*dt/rk_coef(s, 4)) + if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(1)%vf) if (model_eqns == model_eqns_6eq .and. (.not. relax)) then @@ -704,6 +708,33 @@ contains end subroutine s_apply_bodyforces + subroutine s_apply_synthetic_turbulence_force(q_cons_vf, q_prim_vf_in, rhs_vf_in, ldt) + + type(scalar_field), dimension(1:sys_size), intent(inout) :: q_cons_vf + type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf_in + type(scalar_field), dimension(1:sys_size), intent(inout) :: rhs_vf_in + real(wp), intent(in) :: ldt !< local dt + integer :: i, j, k, l + + call nvtxStartRange("RHS-SYNTHETICFORCE") + call s_compute_synthetic_forces_rhs(q_prim_vf_in, q_cons_vf, rhs_vf_in) + + $:GPU_PARALLEL_LOOP(collapse=4) + do i = eqn_idx%mom%beg, eqn_idx%E + do l = 0, p + do k = 0, n + do j = 0, m + q_cons_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) + ldt*rhs_vf_in(i)%sf(j, k, l) + end do + end do + end do + end do + $:END_GPU_PARALLEL_LOOP() + + call nvtxEndRange + + end subroutine s_apply_synthetic_turbulence_force + !> Update immersed boundary positions and velocities at the current Runge-Kutta stage subroutine s_propagate_immersed_boundaries(s) diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 48ce9d74e7..8a04f512cd 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -787,6 +787,23 @@ def _load(): _r(f"p_{d}", REAL, math=r"\f$\phi_" + d + r"\f$") _r(f"bf_{d}", LOG) + # Synthetic turbulence forcing + _r("synthetic_turbulence", LOG, {"synthetic_turbulence"}) + _r("synth_seed", INT, {"synthetic_turbulence"}) + _r("synth_n_shells", INT, {"synthetic_turbulence"}) + _r("num_turbulent_sources", INT, {"synthetic_turbulence"}) + _r("synth_U_inf", REAL, {"synthetic_turbulence"}, math=r"\f$U_\infty\f$") + NSS = _fc("num_synth_shells_max", 50) + NTS = _fc("num_turb_sources_max", 10) + for i in range(1, NSS + 1): + _r(f"synth_n_waves_per_shell({i})", INT, {"synthetic_turbulence"}) + _r(f"synth_k_shell({i})", REAL, {"synthetic_turbulence"}, math=r"\f$k_s\f$") + _r(f"synth_amp_shell({i})", REAL, {"synthetic_turbulence"}, math=r"\f$A_s\f$") + for i in range(1, NTS + 1): + for d in range(1, 4): + _r(f"turb_pos({i},{d})", REAL, {"synthetic_turbulence"}) + _r(f"synth_L({i},{d})", REAL, {"synthetic_turbulence"}) + # INDEXED PARAMETERS # patch_icpp (10 patches) @@ -1270,6 +1287,16 @@ def _nv(targets: set, *names: str) -> None: "g_x", "g_y", "g_z", + "synthetic_turbulence", + "synth_seed", + "synth_n_shells", + "num_turbulent_sources", + "synth_U_inf", + "synth_n_waves_per_shell", + "synth_k_shell", + "synth_amp_shell", + "turb_pos", + "synth_L", "collision_model", "coefficient_of_restitution", "collision_time", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index f9e7e13863..cfe9d14876 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -232,6 +232,12 @@ "p_x": "Body force phase in x-direction", "p_y": "Body force phase in y-direction", "p_z": "Body force phase in z-direction", + # Synthetic turbulence + "synthetic_turbulence": "Enable synthetic turbulence forcing", + "synth_seed": "Random seed for wave vector generation", + "synth_n_shells": "Number of energy shells for synthetic turbulence", + "num_turbulent_sources": "Number of Gaussian forcing zones", + "synth_U_inf": "Advection velocity for synthetic turbulence field", # Output flags "mom_wrt": "Write momentum to database", "flux_wrt": "Write flux data",