Fix magic number and case file checker

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) |
@@ -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.
 
 - `%%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.
 

src/common/m_constants.fpp

@@ -23,7 +23,14 @@ module m_constants
     integer, parameter :: pathlen_max = 400
     integer, parameter :: nnode = 4    !< Number of QBMM nodes
     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 :: acoustic_spatial_support_width = 20d0 !< 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 2004 )
+    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,7 +136,7 @@ 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)), dimension(1:100) :: f_BB, period_BB, sl_BB, bwid_BB, ffre_BB, phi_rn
+        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
@@ -178,13 +178,11 @@ contains
             sum_BB = 0d0
 
             !$acc loop
-            do k = 1, 100
-                f_BB(k) = 500d0 + k*100d0 ! Discrete frequency specturm center
-                period_BB(k) = 1d0/f_BB(k)
-                sl_BB(k) = 20d0*mag(ai) + k*mag(ai)/10 ! Spectral level at each frequency
-                bwid_BB(k) = 100d0 ! Bandwidth
-                ffre_BB(k) = dsqrt((2d0*sl_BB(k)*bwid_BB(k)))*cos((sim_time)*2.d0*pi/period_BB(k) + 2d0*pi*phi_rn(k))
-                sum_BB = sum_BB + ffre_BB(k)
+            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)
@@ -367,10 +365,8 @@ contains
                 source = mag(ai)*sine_wave*1d2
             end if
 
-        elseif (pulse(ai) == 4) then
-            ! TO DO: delay broadband acoustic source
+        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 @@ contains
                 "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)), &