Add Broadband Acoustic Src

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) |
@@ -576,7 +576,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.
 

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,92 @@
+#!/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'                  : 50000,
+'t_step_save'                  : 500,
+# ==========================================================================
+
+# 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'              : 10.,
+# ==========================================================================
+
+# Fluids Physical Parameters ===============================================
+'fluid_pp(1)%gamma'            : 1.E+00/(4.4E+00-1.E+00),
+'fluid_pp(1)%pi_inf'           : 0,
+# ==========================================================================
+}))
+
+# ==============================================================================

src/common/m_constants.fpp

@@ -26,4 +26,11 @@ module m_constants
     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
 
+    !< Broadband acoustic source constant (reference: Tam et al. JSV 2005)
+    integer, parameter :: num_broadband_freq = 100 !< The number of sine wave frequencies in broadband acoustic source.
+    real(kind(0d0)), parameter :: broadband_freq_lowest = 500d0 !< The lower sine wave frequency bound in broadband acoustic source.
+    real(kind(0d0)), parameter :: broadband_bandwidth = 100d0 !< The bandwidth of the discretized sine wave frequencies.
+    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/simulation/m_acoustic_src.fpp

@@ -136,6 +136,8 @@
         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)), dimension(1:num_broadband_freq) :: period_BB, sl_BB, ffre_BB, phi_rn ! broadband source variables
+        real(kind(0d0)) :: sum_BB
 
         integer :: i, j, k, l, q !< generic loop variables
         integer :: ai !< acoustic source index
@@ -171,6 +173,18 @@
 
             num_points = source_spatials_num_points(ai) ! Use scalar to force firstprivate to prevent GPU bug
 
+            call random_number(phi_rn(1:100))
+            call s_mpi_send_random_number(phi_rn)
+            sum_BB = 0d0
+
+            !$acc loop
+            do k = 1, num_broadband_freq
+                period_BB(k) = 1d0/(broadband_freq_lowest + k*broadband_bandwidth) ! Acoustic period of the wave at each discrete frequency
+                sl_BB(k) = broadband_spectral_level_constant*mag(ai) + k*mag(ai)/broadband_spectral_level_growth_rate ! Spectral level at each frequency
+                ffre_BB(k) = dsqrt((2d0*sl_BB(k)*broadband_bandwidth))*cos((sim_time)*2.d0*pi/period_BB(k) + 2d0*pi*phi_rn(k)) ! Source term corresponding to each frequencies
+                sum_BB = sum_BB + ffre_BB(k) ! Total source term for the broadband wave
+            end do
+
             !$acc parallel loop gang vector default(present) private(myalpha, myalpha_rho)
             do i = 1, num_points
                 j = source_spatials(ai)%coord(1, i)
@@ -220,7 +234,7 @@
                 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 +271,7 @@
                     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 +311,10 @@
     !! @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 +365,8 @@
                 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
 

src/simulation/m_checker.fpp

@@ -164,8 +164,8 @@
                 "acoustic("//trim(jStr)//")%mag must be specified")
             @:PROHIBIT(acoustic(j)%pulse == dflt_int, &
                 "acoustic("//trim(jStr)//")%pulse must be specified")
-            @:PROHIBIT(.not. any(acoustic(j)%pulse == (/1, 2, 3/)), &
-                "Only acoustic("//trim(jStr)//")%pulse = 1, 2, or 3 is allowed")
+            @:PROHIBIT(.not. any(acoustic(j)%pulse == (/1, 2, 3, 4/)), &
+                "Only acoustic("//trim(jStr)//")%pulse = 1, 2, 3 or 4 is allowed")
 
             @:PROHIBIT(any(acoustic(j)%pulse == (/1, 3/)) .and. &
                 (f_is_default(acoustic(j)%frequency) .eqv. f_is_default(acoustic(j)%wavelength)), &

src/simulation/m_mpi_proxy.fpp

@@ -2319,6 +2319,13 @@ contains
 
     end subroutine s_mpi_sendrecv_capilary_variables_buffers
 
+    subroutine s_mpi_send_random_number(phi_rn)
+        real(kind(0d0)), dimension(1:100) :: phi_rn
+#ifdef MFC_MPI
+        call MPI_BCAST(phi_rn, size(phi_rn), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+#endif
+    end subroutine s_mpi_send_random_number
+
     !> Module deallocation and/or disassociation procedures
     subroutine s_finalize_mpi_proxy_module