Add Broadband Acoustic Src (#706)

Co-authored-by: Spencer Bryngelson <[email protected]>

Author: haochey

docs/documentation/case.md

@@ -545,7 +545,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 | `acoustic(i)%%support`                | Integer | Geometry of spatial support for the acoustic source |
 | `acoustic(i)%%dipole`                 | Logical | Dipole source activation (optional; default = false -> monopole) |
 | `acoustic(i)%%loc(j)`                 | Real    | $j$-th coordinate of the point that defines the acoustic source location |
-| `acoustic(i)%%pulse`                  | Integer | Acoustic wave form: [1] Sine [2] Gaussian [3] Square |
+| `acoustic(i)%%pulse`                  | Integer | Acoustic wave form: [1] Sine [2] Gaussian [3] Square [4] Broadband  |
 | `acoustic(i)%%npulse`                 | Real    | Number of pulse cycles |
 | `acoustic(i)%%mag`                    | Real    | Pulse magnitude	|
 | `acoustic(i)%%frequency`              | Real    | Sine/Square - Frequency of the acoustic wave  (exclusive) |
@@ -563,6 +563,9 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 | `acoustic(i)%%element_spacing_angle`  | Real    | 2D Transducer array - Spacing angle (in rad) between adjacent transducer elements |
 | `acoustic(i)%%element_polygon_ratio`  | Real    | 3D Transducer array - Ratio of polygon side length to transducer element radius |
 | `acoustic(i)%%rotate_angle`           | Real    | 3D Transducer array - Rotation angle of the transducer array (optional; default = 0) |
+| `acoustic(i)%%bb_num_freq`            | integer | Number of frequencies in broadband wave |
+| `acoustic(i)%%bb_bandwidth`           | Real    | The bandwidth of each frequency in the broadband wave |
+| `acoustic(i)%%bb_lowest_freq`         | Real    | The lower frequency bound of the broadband wave |
 
 Details of the transducer acoustic source model can be found in [Maeda and Colonius (2017)](references.md#Maeda17).
 
@@ -576,7 +579,7 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon
 
 - `%%loc(j)` specifies the location of the acoustic source in the $j$-th coordinate direction. For planer support, the location defines midpoint of the source plane. For transducer arrays, the location defines the center of the transducer or transducer array (not the focal point; for 3D it's the tip of the spherical cap, for 2D it's the tip of the arc). 
 
-- `%%pulse` specifies the acoustic wave form. `%%pulse = 1`, `2`, and `3` correspond to sinusoidal wave, Gaussian wave, and square wave, respectively.
+- `%%pulse` specifies the acoustic wave form. `%%pulse = 1`, `2`, `3` and `4` correspond to sinusoidal wave, Gaussian wave, square wave and broadband wave, respectively. The implementation of the broadband wave is based on [Tam (2005)](references.md#Tam05)
 
 - `%%npulse` specifies the number of cycles of the acoustic wave generated. Only applies to `%%pulse = 1 and 3` (sine and square waves), and must be an integer for non-planar waves.
 
@@ -608,6 +611,12 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon
 
 - `%%rotate_angle` specifies the rotation angle of the 3D circular transducer array along the x-axis (principal axis). It is optional and defaults to 0.
 
+- `%%bb_num_freq` specifies the number discretized frequencies in the broadband acoustic wave. If `%%bb_num_freq` is 1, the acoustic wave will be a discrete tone (i.e. single frequency sine wave).
+
+- `%%bb_bandwidth` specifies the bandwidth of the discretized frequencies.
+
+- `%%bb_lowest_freq` specifies the lower frequency bound of the broadband acoustic wave. The upper frequency bound will be calculated as `%%bb_lowest_freq + %%bb_num_freq * %%bb_bandwidth`. The wave is no longer broadband below the lower bound and above the upper bound.
+
 ### 9. Ensemble-Averaged Bubble Model
 
 | Parameter      | Type    | Description                                    |

docs/documentation/references.md

@@ -40,6 +40,8 @@
 
 - <a id="Suresh97">Suresh, A. and Huynh, H. (1997). Accurate monotonicity-preserving schemes with runge–kutta time stepping. Journal of Computational Physics, 136(1):83–99.</a>
 
+- <a id="Tam05">Tam, C. K., Ju, H., Jones, M. G., Watson, W. R., and Parrott, T. L. (2005). A computational and experimental study of slit resonators. Journal of Sound and Vibration, 284(3-5), 947-984.</a>
+
 - <a id="Thompson87">Thompson, K. W. (1987). Time dependent boundary conditions for hyperbolic systems. Journal of computational physics, 68(1):1–24.</a>
 
 - <a id="Thompson90">Thompson, K. W. (1990). Time-dependent boundary conditions for hyperbolic systems, ii. Journal of computational physics, 89(2):439–461.</a>

examples/2D_acoustic_broadband/case.py

@@ -0,0 +1,95 @@
+#!/usr/bin/env python3
+
+import json, math
+
+print(json.dumps({
+# Logistics ================================================================
+'run_time_info' : 'T',
+# ==========================================================================
+
+# Computational Domain Parameters ==========================================
+'x_domain%beg'                 : 0,
+'x_domain%end'                 : 0.3,
+'y_domain%beg'                 : 0,
+'y_domain%end'                 : 0.1,
+'m'                            : 299,
+'n'                            : 99,
+'p'                            : 0,
+'dt'                           : 5e-7,
+'t_step_start'                 : 0,
+'t_step_stop'                  : 1000,
+'t_step_save'                  : 10,
+# ==========================================================================
+
+# Simulation Algorithm Parameters ==========================================
+'num_patches'                  : 1,
+'model_eqns'                   : 2,
+'alt_soundspeed'               : 'F',
+'num_fluids'                   : 1,
+'mpp_lim'                      : 'F',
+'mixture_err'                  : 'F',
+'time_stepper'                 : 3,
+'weno_order'                   : 5,
+'weno_eps'                     : 1.E-16,
+'teno'                         : 'T',
+'teno_CT'                      : 1E-8,
+'null_weights'                 : 'F',
+'mp_weno'                      : 'F',
+'riemann_solver'               : 2,
+'wave_speeds'                  : 1,
+'avg_state'                    : 2,
+'bc_x%beg'                     : -6,
+'bc_x%end'                     : -6,
+'bc_y%beg'                     : -6,
+'bc_y%end'                     : -6,
+# ==========================================================================
+
+# Formatted Database Files Structure Parameters ============================
+'format'                       : 1,
+'precision'                    : 2,
+'prim_vars_wrt'                :'T',
+'parallel_io'                  :'T',
+'probe_wrt'                    :'T',
+'fd_order'                     : 2,
+'num_probes'                   : 1,
+'probe(1)%x'                   : 0.13,
+'probe(1)%y'                   : 0.05,
+
+# ==========================================================================
+
+# Patch 1 Liquid ===========================================================
+'patch_icpp(1)%geometry'       : 3,
+'patch_icpp(1)%x_centroid'     : 0.15,
+'patch_icpp(1)%y_centroid'     : 0.05,
+'patch_icpp(1)%length_x'       : 0.3,
+'patch_icpp(1)%length_y'       : 0.1,
+'patch_icpp(1)%vel(1)'         : 0.0,
+'patch_icpp(1)%vel(2)'         : 0.0,
+'patch_icpp(1)%pres'           : 1E+05,
+'patch_icpp(1)%alpha_rho(1)'   : 1.19,
+'patch_icpp(1)%alpha(1)'       : 1.0,
+# ==========================================================================
+
+# Acoustic source ==========================================================
+'acoustic_source'              : 'T',
+'num_source'                   : 1,
+'acoustic(1)%support'          : 2,
+'acoustic(1)%loc(1)'           : 0.1,
+'acoustic(1)%loc(2)'           : 0.05,
+'acoustic(1)%dir'              : math.pi * 0,
+'acoustic(1)%length'           : 0.1,
+'acoustic(1)%pulse'            : 4,
+'acoustic(1)%npulse'           : 1000000,
+'acoustic(1)%mag'              : 1000.,
+'acoustic(1)%bb_num_freq'      : 100,
+'acoustic(1)%bb_lowest_freq'   : 500.,
+'acoustic(1)%bb_bandwidth'     : 500.,
+# ==========================================================================
+
+# Fluids Physical Parameters ===============================================
+'fluid_pp(1)%gamma'            : 1.E+00/(1.4E+00-1.E+00),
+'fluid_pp(1)%pi_inf'           : 0,
+# ==========================================================================
+}))
+
+# ==============================================================================

src/common/m_constants.fpp

@@ -25,5 +25,7 @@ module m_constants
     real(kind(0d0)), parameter :: capillary_cutoff = 1e-6 !< color function gradient magnitude at which to apply the surface tension fluxes
     real(kind(0d0)), parameter :: acoustic_spatial_support_width = 2.5d0 !< Spatial support width of acoustic source, used in s_source_spatial
     real(kind(0d0)), parameter :: dflt_vcfl_dt = 100d0 !< value of vcfl_dt when viscosity is off for computing adaptive timestep size
+    real(kind(0d0)), parameter :: broadband_spectral_level_constant = 20d0 !< The constant to scale the spectral level at the lower frequency bound
+    real(kind(0d0)), parameter :: broadband_spectral_level_growth_rate = 10d0 !< The spectral level constant to correct the magnitude at each frqeuency to ensure the source is overall broadband
 
 end module m_constants

src/common/m_derived_types.fpp

@@ -298,8 +298,11 @@ module m_derived_types
         real(kind(0d0)) :: element_spacing_angle !< Spacing between aperture elements in 2D acoustic array
         real(kind(0d0)) :: element_polygon_ratio !< Ratio of aperture element diameter to side length of polygon connecting their centers, in 3D acoustic array
         real(kind(0d0)) :: rotate_angle !< Angle of rotation of the entire circular 3D acoustic array
+        real(kind(0d0)) :: bb_bandwidth !< Bandwidth of each frequency in broadband wave
+        real(kind(0d0)) :: bb_lowest_freq !< The lower frequency bound of broadband wave
         integer :: num_elements !< Number of elements in the acoustic array
         integer :: element_on !< Element in the acoustic array to turn on
+        integer :: bb_num_freq !< Number of frequencies in the broadband wave
     end type acoustic_parameters
 
     !> Acoustic source source_spatial pre-calculated values

src/simulation/m_acoustic_src.fpp

@@ -41,8 +41,11 @@ module m_acoustic_src
     real(kind(0d0)), allocatable, dimension(:) :: element_spacing_angle, element_polygon_ratio, rotate_angle
     !$acc declare create(element_spacing_angle, element_polygon_ratio, rotate_angle)
 
-    integer, allocatable, dimension(:) :: num_elements, element_on
-    !$acc declare create(num_elements, element_on)
+    real(kind(0d0)), allocatable, dimension(:) :: bb_bandwidth, bb_lowest_freq
+    !$acc declare create(bb_bandwidth, bb_lowest_freq)
+
+    integer, allocatable, dimension(:) :: num_elements, element_on, bb_num_freq
+    !$acc declare create(num_elements, element_on, bb_num_freq)
 
     !> @name Acoustic source terms
     !> @{
@@ -63,7 +66,8 @@ contains
     subroutine s_initialize_acoustic_src
         integer :: i, j !< generic loop variables
 
-        @:ALLOCATE(loc_acoustic(1:3, 1:num_source), mag(1:num_source), dipole(1:num_source), support(1:num_source), length(1:num_source), height(1:num_source), wavelength(1:num_source), frequency(1:num_source), gauss_sigma_dist(1:num_source), gauss_sigma_time(1:num_source), foc_length(1:num_source), aperture(1:num_source), npulse(1:num_source), pulse(1:num_source), dir(1:num_source), delay(1:num_source), element_polygon_ratio(1:num_source), rotate_angle(1:num_source), element_spacing_angle(1:num_source), num_elements(1:num_source), element_on(1:num_source))
+        @:ALLOCATE(loc_acoustic(1:3, 1:num_source), mag(1:num_source), dipole(1:num_source), support(1:num_source), length(1:num_source), height(1:num_source), wavelength(1:num_source), frequency(1:num_source), gauss_sigma_dist(1:num_source), gauss_sigma_time(1:num_source), foc_length(1:num_source), aperture(1:num_source), npulse(1:num_source), pulse(1:num_source), dir(1:num_source), delay(1:num_source), element_polygon_ratio(1:num_source), rotate_angle(1:num_source), element_spacing_angle(1:num_source), num_elements(1:num_source), element_on(1:num_source), bb_num_freq(1:num_source), bb_bandwidth(1:num_source), bb_lowest_freq(1:num_source))
+
         do i = 1, num_source
             do j = 1, 3
                 loc_acoustic(j, i) = acoustic(i)%loc(j)
@@ -85,6 +89,10 @@ contains
             element_spacing_angle(i) = acoustic(i)%element_spacing_angle
             element_polygon_ratio(i) = acoustic(i)%element_polygon_ratio
             num_elements(i) = acoustic(i)%num_elements
+            bb_num_freq(i) = acoustic(i)%bb_num_freq
+            bb_bandwidth(i) = acoustic(i)%bb_bandwidth
+            bb_lowest_freq(i) = acoustic(i)%bb_lowest_freq
+
             if (acoustic(i)%element_on == dflt_int) then
                 element_on(i) = 0
             else
@@ -101,7 +109,7 @@ contains
                 delay(i) = acoustic(i)%delay
             end if
         end do
-        !$acc update device(loc_acoustic, mag, dipole, support, length, height, wavelength, frequency, gauss_sigma_dist, gauss_sigma_time, foc_length, aperture, npulse, pulse, dir, delay, element_polygon_ratio, rotate_angle, element_spacing_angle, num_elements, element_on)
+        !$acc update device(loc_acoustic, mag, dipole, support, length, height, wavelength, frequency, gauss_sigma_dist, gauss_sigma_time, foc_length, aperture, npulse, pulse, dir, delay, element_polygon_ratio, rotate_angle, element_spacing_angle, num_elements, element_on, bb_num_freq, bb_bandwidth, bb_lowest_freq)
 
         @:ALLOCATE(mass_src(0:m, 0:n, 0:p))
         @:ALLOCATE(mom_src(1:num_dims, 0:m, 0:n, 0:p))
@@ -136,6 +144,11 @@ contains
         real(kind(0d0)) :: frequency_local, gauss_sigma_time_local
         real(kind(0d0)) :: mass_src_diff, mom_src_diff
         real(kind(0d0)) :: source_temporal
+        real(kind(0d0)) :: period_BB !< period of each sine wave in broadband source
+        real(kind(0d0)) :: sl_BB !< spectral level at each frequency
+        real(kind(0d0)) :: ffre_BB !< source term corresponding to each frequency
+        real(kind(0d0)) :: sum_BB !< total source term for the broadband wave
+        real(kind(0d0)), allocatable, dimension(:) :: phi_rn !< random phase shift for each frequency
 
         integer :: i, j, k, l, q !< generic loop variables
         integer :: ai !< acoustic source index
@@ -171,6 +184,35 @@ contains
 
             num_points = source_spatials_num_points(ai) ! Use scalar to force firstprivate to prevent GPU bug
 
+            ! Calculate the broadband source
+            period_BB = 0d0
+            sl_BB = 0d0
+            ffre_BB = 0d0
+            sum_BB = 0d0
+
+            ! Allocate buffers for random phase shift
+            allocate (phi_rn(1:bb_num_freq(ai)))
+
+            if (pulse(ai) == 4) then
+                call random_number(phi_rn(1:bb_num_freq(ai)))
+                ! Ensure all the ranks have the same random phase shift
+                call s_mpi_send_random_number(phi_rn, bb_num_freq(ai))
+            end if
+
+            !$acc loop reduction(+:sum_BB)
+            do k = 1, bb_num_freq(ai)
+                ! Acoustic period of the wave at each discrete frequency
+                period_BB = 1d0/(bb_lowest_freq(ai) + k*bb_bandwidth(ai))
+                ! Spectral level at each frequency
+                sl_BB = broadband_spectral_level_constant*mag(ai) + k*mag(ai)/broadband_spectral_level_growth_rate
+                ! Source term corresponding to each frequencies
+                ffre_BB = dsqrt((2d0*sl_BB*bb_bandwidth(ai)))*cos((sim_time)*2d0*pi/period_BB + 2d0*pi*phi_rn(k))
+                ! Sum up the source term of each frequency to obtain the total source term for broadband wave
+                sum_BB = sum_BB + ffre_BB
+            end do
+
+            deallocate (phi_rn)
+
             !$acc parallel loop gang vector default(present) private(myalpha, myalpha_rho)
             do i = 1, num_points
                 j = source_spatials(ai)%coord(1, i)
@@ -190,7 +232,7 @@ contains
 
                 if (bubbles) then
                     if (num_fluids > 2) then
-                        !$acc loop reduction(+:myRho,B_tait,small_gamma)
+                        !$acc loop seq
                         do q = 1, num_fluids - 1
                             myRho = myRho + myalpha_rho(q)
                             B_tait = B_tait + myalpha(q)*pi_infs(q)
@@ -204,7 +246,7 @@ contains
                 end if
 
                 if ((.not. bubbles) .or. (mpp_lim .and. (num_fluids > 2))) then
-                    !$acc loop reduction(+:myRho,B_tait,small_gamma)
+                    !$acc loop seq
                     do q = 1, num_fluids
                         myRho = myRho + myalpha_rho(q)
                         B_tait = B_tait + myalpha(q)*pi_infs(q)
@@ -220,7 +262,7 @@ contains
                 if (pulse(ai) == 2) gauss_sigma_time_local = f_gauss_sigma_time_local(gauss_conv_flag, ai, c)
 
                 ! Update momentum source term
-                call s_source_temporal(sim_time, c, ai, mom_label, frequency_local, gauss_sigma_time_local, source_temporal)
+                call s_source_temporal(sim_time, c, ai, mom_label, frequency_local, gauss_sigma_time_local, source_temporal, sum_BB)
                 mom_src_diff = source_temporal*source_spatials(ai)%val(i)
 
                 if (dipole(ai)) then ! Double amplitude & No momentum source term (only works for Planar)
@@ -257,7 +299,7 @@ contains
                     mass_src_diff = mom_src_diff/c
                 else ! Spherical or cylindrical support
                     ! Mass source term must be calculated differently using a correction term for spherical and cylindrical support
-                    call s_source_temporal(sim_time, c, ai, mass_label, frequency_local, gauss_sigma_time_local, source_temporal)
+                    call s_source_temporal(sim_time, c, ai, mass_label, frequency_local, gauss_sigma_time_local, source_temporal, sum_BB)
                     mass_src_diff = source_temporal*source_spatials(ai)%val(i)
                 end if
                 mass_src(j, k, l) = mass_src(j, k, l) + mass_src_diff
@@ -297,10 +339,10 @@ contains
     !! @param frequency_local Frequency at the spatial location for sine and square waves
     !! @param gauss_sigma_time_local sigma in time for Gaussian pulse
     !! @param source Source term amplitude
-    subroutine s_source_temporal(sim_time, c, ai, term_index, frequency_local, gauss_sigma_time_local, source)
+    subroutine s_source_temporal(sim_time, c, ai, term_index, frequency_local, gauss_sigma_time_local, source, sum_BB)
         !$acc routine seq
         integer, intent(in) :: ai, term_index
-        real(kind(0d0)), intent(in) :: sim_time, c
+        real(kind(0d0)), intent(in) :: sim_time, c, sum_BB
         real(kind(0d0)), intent(in) :: frequency_local, gauss_sigma_time_local
         real(kind(0d0)), intent(out) :: source
 
@@ -351,6 +393,8 @@ contains
                 source = mag(ai)*sine_wave*1d2
             end if
 
+        elseif (pulse(ai) == 4) then ! Broadband wave
+            source = sum_BB
         end if
     end subroutine s_source_temporal