Add broadband frequency input parameters

Author: haochey

docs/documentation/case.md

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

examples/2D_acoustic_broadband/case.py

@@ -17,8 +17,8 @@
 'p'                            : 0,
 'dt'                           : 5e-7,
 't_step_start'                 : 0,
-'t_step_stop'                  : 50000,
-'t_step_save'                  : 500,
+'t_step_stop'                  : 1000,
+'t_step_save'                  : 10,
 # ==========================================================================
 
 # Simulation Algorithm Parameters ==========================================
@@ -80,11 +80,14 @@
 'acoustic(1)%length'           : 0.1,
 'acoustic(1)%pulse'            : 4,
 'acoustic(1)%npulse'           : 1000000,
-'acoustic(1)%mag'              : 10.,
+'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/(4.4E+00-1.E+00),
+'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,11 +25,6 @@ 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
-
-    !< 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
 

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,7 @@ contains
     subroutine s_initialize_acoustic_src
         integer :: i, j !< generic loop variables
 
-        @:ALLOCATE_GLOBAL(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_GLOBAL(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 +88,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 +108,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_GLOBAL(mass_src(0:m, 0:n, 0:p))
         @:ALLOCATE_GLOBAL(mom_src(1:num_dims, 0:m, 0:n, 0:p))
@@ -136,11 +143,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)), dimension(1:num_broadband_freq) :: period_BB !< period of each sine wave in broadband source
-        real(kind(0d0)), dimension(1:num_broadband_freq) :: sl_BB !< spectral level at each frequency
-        real(kind(0d0)), dimension(1:num_broadband_freq) :: ffre_BB !< source term corresponding to each frequency
-        real(kind(0d0)), dimension(1:num_broadband_freq) :: phi_rn !< random phase shift for each frequency
+        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
@@ -176,23 +183,33 @@ contains
 
             num_points = source_spatials_num_points(ai) ! Use scalar to force firstprivate to prevent GPU bug
 
-            call random_number(phi_rn(1:num_broadband_freq))
-            ! Ensure all the ranks have the same random phase shift
-            call s_mpi_send_random_number(phi_rn)
+            ! 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)))
+
+            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))
+
             !$acc loop reduction(+:sum_BB)
-            do k = 1, num_broadband_freq
+            do k = 1, bb_num_freq(ai)
                 ! Acoustic period of the wave at each discrete frequency
-                period_BB(k) = 1d0/(broadband_freq_lowest + k*broadband_bandwidth)
+                period_BB = 1d0/(bb_lowest_freq(ai) + k*bb_bandwidth(ai))
                 ! Spectral level at each frequency
-                sl_BB(k) = broadband_spectral_level_constant*mag(ai) + k*mag(ai)/broadband_spectral_level_growth_rate
+                sl_BB = broadband_spectral_level_constant*mag(ai) + k*mag(ai)/broadband_spectral_level_growth_rate
                 ! Source term corresponding to each frequencies
-                ffre_BB(k) = dsqrt((2d0*sl_BB(k)*broadband_bandwidth))*cos((sim_time)*2d0*pi/period_BB(k) + 2d0*pi*phi_rn(k))
+                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(k)
+                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)

src/simulation/m_checker.fpp

@@ -175,6 +175,12 @@ contains
                 (f_is_default(acoustic(j)%gauss_sigma_time) .eqv. f_is_default(acoustic(j)%gauss_sigma_dist)), &
                 "One and only one of acoustic("//trim(jStr)//")%gauss_sigma_time "// &
                 "or acoustic("//trim(jStr)//")%gauss_sigma_dist must be specified for pulse = 2")
+            @:PROHIBIT(acoustic(j)%pulse == 4 .and. acoustic(j)%bb_num_freq == dflt_int, &
+                "The number of broadband frequencies acoustic("//trim(jStr)//")%bb_num_freq must be specified for pulse = 4")
+            @:PROHIBIT(acoustic(j)%pulse == 4 .and. f_is_default(acoustic(j)%bb_bandwidth), &
+                "The broadband wave band width acoustic("//trim(jStr)//")%bb_bandwidth must be specified for pulse = 4")
+            @:PROHIBIT(acoustic(j)%pulse == 4 .and. f_is_default(acoustic(j)%bb_lowest_freq), &
+                "The broadband wave lower frequency bound acoustic("//trim(jStr)//")%bb_lowest_freq must be specified for pulse = 4")
 
             @:PROHIBIT(f_is_default(acoustic(j)%npulse), &
                 "acoustic("//trim(jStr)//")%npulse must be specified")

src/simulation/m_global_parameters.fpp

@@ -644,6 +644,12 @@ contains
             acoustic(j)%rotate_angle = dflt_real
             acoustic(j)%num_elements = dflt_int
             acoustic(j)%element_on = dflt_int
+            acoustic(j)%rotate_angle = dflt_real
+            acoustic(j)%num_elements = dflt_int
+            acoustic(j)%element_on = dflt_int
+            acoustic(j)%bb_num_freq = dflt_int
+            acoustic(j)%bb_lowest_freq = dflt_real
+            acoustic(j)%bb_bandwidth = dflt_real
         end do
 
         fd_order = dflt_int

src/simulation/m_mpi_proxy.fpp

@@ -252,14 +252,15 @@ contains
 
             call MPI_BCAST(acoustic(j)%dipole, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
 
-            #:for VAR in [ 'pulse', 'support', 'num_elements', 'element_on' ]
+            #:for VAR in [ 'pulse', 'support', 'num_elements', 'element_on', 'bb_num_freq' ]
                 call MPI_BCAST(acoustic(j)%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 
             #:for VAR in [ 'mag', 'length', 'height', &
                 'wavelength', 'frequency', 'gauss_sigma_dist', 'gauss_sigma_time', &
                 'npulse', 'dir', 'delay', 'foc_length', 'aperture', &
-                'element_spacing_angle', 'element_polygon_ratio', 'rotate_angle' ]
+                'element_spacing_angle', 'element_polygon_ratio', 'rotate_angle', &
+                'bb_bandwidth', 'bb_lowest_freq' ]
                 call MPI_BCAST(acoustic(j)%${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 
@@ -2332,10 +2333,11 @@ contains
 
     end subroutine s_mpi_sendrecv_capilary_variables_buffers
 
-    subroutine s_mpi_send_random_number(phi_rn)
-        real(kind(0d0)), intent(inout), dimension(1:num_broadband_freq) :: phi_rn
+    subroutine s_mpi_send_random_number(phi_rn, num_freq)
+        integer, intent(in) :: num_freq
+        real(kind(0d0)), intent(inout), dimension(1:num_freq) :: phi_rn
 #ifdef MFC_MPI
-        call MPI_BCAST(phi_rn, size(phi_rn), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        call MPI_BCAST(phi_rn, num_freq, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
 #endif
     end subroutine s_mpi_send_random_number
 

toolchain/mfc/run/case_dicts.py

@@ -279,7 +279,8 @@ def analytic(self):
         SIMULATION[f"fluid_pp({f_id})%Re({re_id})"] = ParamType.REAL
 
     for mono_id in range(1,4+1):
-        for int_attr in ["pulse", "support", "num_elements", "element_on"]:
+        for int_attr in ["pulse", "support", "num_elements", "element_on", 
+                         "bb_num_freq"]:
             SIMULATION[f"acoustic({mono_id})%{int_attr}"] = ParamType.INT
 
         SIMULATION[f"acoustic({mono_id})%dipole"] = ParamType.LOG
@@ -288,7 +289,7 @@ def analytic(self):
                           "gauss_sigma_dist", "gauss_sigma_time", "npulse",
                           "dir", "delay", "foc_length", "aperture",
                           "element_spacing_angle", "element_polygon_ratio",
-                          "rotate_angle"]:
+                          "rotate_angle", "bb_bandwidth", "bb_lowest_freq"]:
             SIMULATION[f"acoustic({mono_id})%{real_attr}"] = ParamType.REAL
 
         for cmp_id in range(1,3+1):