pass all GPU tests locally

docs/documentation/case.md

@@ -815,7 +815,7 @@ For all acoustic support types:
 - `%%npulse == 2` requires `%%delay` to be set.
 
 Description of the acoustic support types:
-- `%%support = 1` specifies an infinite source plane that is normal to the $x$-axis and intersects with the axis at $x=$ `%%loc(1)` in 1D simulation.
+- `%%support = 1` specifies an infinite source plane that is normal to the $x$-axis and intersects with the axis at $x=$ `%%loc(1)` in 1D simulation. `%%dir > 0` specifies a rightward propagating wave, and `%%dir < 0` specifies a leftward propagating wave. `%%dir = 0` is not allowed.
 
 - `%%support = 2` specifies a semi-infinite source plane in 2D simulation.
 The midplane location is [`%%loc(1)`, `%%loc(2)`] and the normal vector is [$\mathrm{cos}$(`%%dir`), $\mathrm{sin}$(`%%dir`)]. The length of the source plane is `%%length`, and the plane is perpendicular to the direction of wave propagation (defined by `%%dir`).
@@ -829,7 +829,7 @@ The midplane location is [`%%loc(1)`, `%%loc(2)`] and the normal vector is [$\ma
 
 - `%%support = 7` specifies a spherical transducer in 3D simulation. The transducer is centered at [`%%loc(1)`, `%%loc(2)`, `%%loc(3)`] with a focal length of `%%foc_length` and an aperture of `%%aperture`. The center location is not the focal point; it is the tip of the spherical cap (intersection of the cap and the x-axis). The aperture is the diameter of the projection of the spherical cap onto the y-z plane. If a semi-sphere is desired, set the aperture to double the focal length. Again, the equation is exact.
 
-- `%%support = 9` specifies an arcuate transducer array in 2D simulation. The total aperture of the array is `%%aperture`, which is similar to `%%support = 5`. The parameters `%%num_elements` and `%%element_spacing_angle` specify the number of transducer elements and the spacing angle. The spacing angle is the angle of the gap betwen adjacent transducer elements in the array. Because the total aperture is set, each transducer element is smaller if the spacing angle is larger. Physically it represents curved panels. Note that similar to `%%support = 5`, the mass source term correction factor is empirically determined to be -0.85.
+- `%%support = 9` specifies an arcuate transducer array in 2D simulation. The total aperture of the array is `%%aperture`, which is similar to `%%support = 5`. The parameters `%%num_elements` and `%%element_spacing_angle` specify the number of transducer elements and the spacing angle. The spacing angle is the angle of the gap between adjacent transducer elements in the array. Because the total aperture is set, each transducer element is smaller if the spacing angle is larger. Physically it represents curved panels. Note that similar to `%%support = 5`, the mass source term correction factor is empirically determined to be -0.85.
 
 - `%%support = 10` specifies an annular transducer array in 2D axisymmetric simulation. It is identical to `%%support = 9` in terms of simulation parameters. It physically represents the a annulus obtained by revolving the arc in `%%support = 9` around the x-axis.
 

src/simulation/m_checker.fpp

@@ -239,13 +239,13 @@ contains
             if (any(acoustic(j)%pulse == (/1, 3/)) .and. &
                 .not. xor(f_is_default(acoustic(j)%frequency), f_is_default(acoustic(j)%wavelength))) then
                 call s_mpi_abort('One and only one of acoustic(i)frequency '// &
-                                    'or acoustic(i)wavelength must be specified '// &
-                                    'for acoustic(i)pulse = 1 or 3. Exiting ...')
+                                 'or acoustic(i)wavelength must be specified '// &
+                                 'for acoustic(i)pulse = 1 or 3. Exiting ...')
             elseif (acoustic(j)%pulse == 2 .and. &
                     .not. xor(f_is_default(acoustic(j)%gauss_sigma_time), f_is_default(acoustic(j)%gauss_sigma_dist))) then
                 call s_mpi_abort('One and only one of acoustic(i)gauss_sigma_time '// &
-                                    'or acoustic(i)gauss_sigma_dist must be specified '// &
-                                    'for acoustic(i)pulse = 2. Exiting ...')
+                                 'or acoustic(i)gauss_sigma_dist must be specified '// &
+                                 'for acoustic(i)pulse = 2. Exiting ...')
             end if
             if (f_is_default(acoustic(j)%npulse)) then
                 call s_mpi_abort('acoustic('//trim(jStr)//')%npulse must be '// &

src/simulation/m_monopole.fpp

@@ -146,7 +146,7 @@ contains
         real(kind(0d0)) :: frequency_local, gauss_sigma_time_local
         real(kind(0d0)) :: mass_src_diff, mom_src_diff
         real(kind(0d0)) :: angle ! Angle with x-axis for mom source term vector
-        real(kind(0d0)) :: xyz_to_r_ratios(3) ! [xyz]/r for mom source term vector
+        real(kind(0d0)) :: x_to_r_ratio, y_to_r_ratio, z_to_r_ratio ! [xyz]/r for mom source term vector
         real(kind(0d0)) :: source_temporal, source_spatial
 
         integer :: i, j, k, l, q !< generic loop variables
@@ -175,14 +175,12 @@ contains
         do ai = 1, num_source
             ! Skip if the pulse has not started yet for sine and square waves
             if (sim_time < delay(ai) .and. (pulse(ai) == 1 .or. pulse(ai) == 3)) cycle
-            
+
             ! Decide if frequency need to be converted from wavelength
-            frequency_local = frequency(ai)
-            gauss_sigma_time_local = gauss_sigma_time(ai)
             freq_conv_flag = f_is_default(frequency(ai))
             gauss_conv_flag = f_is_default(gauss_sigma_time(ai))
-            
-            !$acc parallel loop collapse(3) gang vector default(present) private(q_cons_local, xyz_to_r_ratios)
+
+            !$acc parallel loop collapse(3) gang vector default(present) private(q_cons_local)
             do l = 0, p
                 do k = 0, n
                     do j = 0, m
@@ -196,15 +194,12 @@ contains
                         call s_compute_speed_of_sound_acoustic_src(q_cons_local, q_prim_local, c, small_gamma)
 
                         ! Wavelength to frequency conversion
-                        if (freq_conv_flag .and. (pulse(ai) == 1 .or. pulse(ai) == 3)) then
-                            frequency_local = c/wavelength(ai)
-                        elseif (gauss_conv_flag .and. pulse(ai) == 2) then
-                            gauss_sigma_time_local = gauss_sigma_dist(ai)/c
-                        end if
+                        if (pulse(ai) == 1 .or. pulse(ai) == 3) frequency_local = f_frequency_local(freq_conv_flag, ai, c)
+                        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_spatial(j, k, l, loc_acoustic(:, ai), ai, source_spatial, angle, xyz_to_r_ratios)
+                        call s_source_spatial(j, k, l, loc_acoustic(:, ai), ai, source_spatial, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio)
                         mom_src_diff = source_temporal*source_spatial
 
                         if (n == 0) then ! 1D
@@ -222,9 +217,9 @@ contains
                                 mom_src(1, j, k, l) = mom_src(1, j, k, l) + mom_src_diff*cos(dir(ai))
                                 mom_src(2, j, k, l) = mom_src(2, j, k, l) + mom_src_diff*sin(dir(ai))
                             else
-                                mom_src(1, j, k, l) = mom_src(1, j, k, l) + mom_src_diff*xyz_to_r_ratios(1)
-                                mom_src(2, j, k, l) = mom_src(2, j, k, l) + mom_src_diff*xyz_to_r_ratios(2)
-                                mom_src(3, j, k, l) = mom_src(3, j, k, l) + mom_src_diff*xyz_to_r_ratios(3)
+                                mom_src(1, j, k, l) = mom_src(1, j, k, l) + mom_src_diff*x_to_r_ratio
+                                mom_src(2, j, k, l) = mom_src(2, j, k, l) + mom_src_diff*y_to_r_ratio
+                                mom_src(3, j, k, l) = mom_src(3, j, k, l) + mom_src_diff*z_to_r_ratio
                             end if
                         end if
 
@@ -233,7 +228,7 @@ contains
                             mass_src_diff = mom_src_diff/c
                         else
                             call s_source_temporal(sim_time, c, ai, mass_label, frequency_local, gauss_sigma_time_local, source_temporal)
-                            call s_source_spatial(j, k, l, loc_acoustic(:, ai), ai, source_spatial, angle, xyz_to_r_ratios)
+                            call s_source_spatial(j, k, l, loc_acoustic(:, ai), ai, source_spatial, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio)
                             mass_src_diff = source_temporal*source_spatial
                         end if
                         mass_src(j, k, l) = mass_src(j, k, l) + mass_src_diff
@@ -333,12 +328,14 @@ contains
     !! @param ai Acoustic source index
     !! @param source Source term amplitude
     !! @param angle Angle of the source term with respect to the x-axis (for 2D or 2D axisymmetric)
-    !! @param xyz_to_r_ratios Ratio of the [xyz]-component of the source term to the magnitude (for 3D)
-    subroutine s_source_spatial(j, k, l, loc, ai, source, angle, xyz_to_r_ratios)
+    !! @param x_to_r_ratio Ratio of the x-component of the source term to the magnitude (for 3D)
+    !! @param y_to_r_ratio Ratio of the y-component of the source term to the magnitude (for 3D) 
+    !! @param z_to_r_ratio Ratio of the z-component of the source term to the magnitude (for 3D)
+    subroutine s_source_spatial(j, k, l, loc, ai, source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio)
         !$acc routine seq
         integer, intent(in) :: j, k, l, ai
         real(kind(0d0)), dimension(3), intent(in) :: loc
-        real(kind(0d0)), intent(out) :: source, angle, xyz_to_r_ratios(3)
+        real(kind(0d0)), intent(out) :: source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio
 
         real(kind(0d0)) :: sig, r(3)
 
@@ -360,9 +357,9 @@ contains
         if (any(support(ai) == (/1, 2, 3, 4/))) then
             call s_source_spatial_planar(ai, sig, r, source)
         elseif (any(support(ai) == (/5, 6, 7/))) then
-            call s_source_spatial_transducer(ai, sig, r, source, angle, xyz_to_r_ratios)
+            call s_source_spatial_transducer(ai, sig, r, source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio)
         elseif (any(support(ai) == (/9, 10, 11/))) then
-            call s_source_spatial_transducer_array(ai, sig, r, source, angle, xyz_to_r_ratios)
+            call s_source_spatial_transducer_array(ai, sig, r, source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio)
         end if
     end subroutine s_source_spatial
 
@@ -399,12 +396,14 @@ contains
     !! @param r Displacement from source to current point
     !! @param source Source term amplitude
     !! @param angle Angle of the source term with respect to the x-axis (for 2D or 2D axisymmetric)
-    !! @param xyz_to_r_ratios Ratio of the [xyz]-component of the source term to the magnitude (for 3D)
-    subroutine s_source_spatial_transducer(ai, sig, r, source, angle, xyz_to_r_ratios)
+    !! @param x_to_r_ratio Ratio of the x-component of the source term to the magnitude (for 3D)
+    !! @param y_to_r_ratio Ratio of the y-component of the source term to the magnitude (for 3D) 
+    !! @param z_to_r_ratio Ratio of the z-component of the source term to the magnitude (for 3D)
+    subroutine s_source_spatial_transducer(ai, sig, r, source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio)
         !$acc routine seq
         integer, intent(in) :: ai
         real(kind(0d0)), intent(in) :: sig, r(3)
-        real(kind(0d0)), intent(out) :: source, angle, xyz_to_r_ratios(3)
+        real(kind(0d0)), intent(out) :: source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio
 
         real(kind(0d0)) :: current_angle, angle_half_aperture, dist, norm
 
@@ -429,9 +428,9 @@ contains
                 source = 1d0/(dsqrt(2d0*pi)*sig/2d0)*dexp(-0.5d0*(dist/(sig/2d0))**2d0)
 
                 norm = dsqrt(r(2)**2d0 + r(3)**2d0 + (foc_length(ai) - r(1))**2d0)
-                xyz_to_r_ratios(1) = -(r(1) - foc_length(ai))/norm
-                xyz_to_r_ratios(2) = -r(2)/norm
-                xyz_to_r_ratios(3) = -r(3)/norm
+                x_to_r_ratio = -(r(1) - foc_length(ai))/norm
+                y_to_r_ratio = -r(2)/norm
+                z_to_r_ratio = -r(3)/norm
             end if
 
         end if
@@ -443,12 +442,14 @@ contains
     !! @param r Displacement from source to current point
     !! @param source Source term amplitude
     !! @param angle Angle of the source term with respect to the x-axis (for 2D or 2D axisymmetric)
-    !! @param xyz_to_r_ratios Ratio of the [xyz]-component of the source term to the magnitude (for 3D)
-    subroutine s_source_spatial_transducer_array(ai, sig, r, source, angle, xyz_to_r_ratios)
+    !! @param x_to_r_ratio Ratio of the x-component of the source term to the magnitude (for 3D)
+    !! @param y_to_r_ratio Ratio of the y-component of the source term to the magnitude (for 3D) 
+    !! @param z_to_r_ratio Ratio of the z-component of the source term to the magnitude (for 3D)
+    subroutine s_source_spatial_transducer_array(ai, sig, r, source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio)
         !$acc routine seq
         integer, intent(in) :: ai
         real(kind(0d0)), intent(in) :: sig, r(3)
-        real(kind(0d0)), intent(out) :: source, angle, xyz_to_r_ratios(3)
+        real(kind(0d0)), intent(out) :: source, angle, x_to_r_ratio, y_to_r_ratio, z_to_r_ratio
 
         integer :: elem, elem_min, elem_max
         real(kind(0d0)) :: current_angle, angle_half_aperture, angle_per_elem, dist
@@ -510,9 +511,9 @@ contains
                     source = dexp(-0.5d0*(dist/(sig/2d0))**2d0)/(dsqrt(2d0*pi)*sig/2d0)
 
                     norm = dsqrt(r(2)**2d0 + r(3)**2d0 + (f - r(1))**2d0)
-                    xyz_to_r_ratios(1) = -(r(1) - f)/norm
-                    xyz_to_r_ratios(2) = -r(2)/norm
-                    xyz_to_r_ratios(3) = -r(3)/norm
+                    x_to_r_ratio = -(r(1) - f)/norm
+                    y_to_r_ratio = -r(2)/norm
+                    z_to_r_ratio = -r(3)/norm
                 end if
 
             end do
@@ -580,4 +581,40 @@ contains
 
     end subroutine s_compute_speed_of_sound_acoustic_src
 
+    !> This function performs wavelength to frequency conversion
+    !! @param freq_conv_flag Determines if frequency is given or wavelength
+    !! @param ai Acoustic source index
+    !! @param c Speed of sound
+    !! @return frequency_local Converted frequency
+    function f_frequency_local(freq_conv_flag, ai, c)
+        logical, intent(in) :: freq_conv_flag
+        integer, intent(in) :: ai
+        real(kind(0d0)), intent(in) :: c
+        real(kind(0d0)) :: f_frequency_local
+    
+        if (freq_conv_flag) then
+            f_frequency_local = c / wavelength(ai)
+        else
+            f_frequency_local = frequency(ai)
+        end if
+    end function f_frequency_local
+
+    !> This function performs Gaussian sigma dist to time conversion
+    !! @param gauss_conv_flag Determines if sigma_dist is given or sigma_time
+    !! @param c Speed of sound
+    !! @param ai Acoustic source index
+    !! @return gauss_sigma_time_local Converted Gaussian sigma time
+    function f_gauss_sigma_time_local(gauss_conv_flag, ai, c)
+        logical, intent(in) :: gauss_conv_flag
+        integer, intent(in) :: ai
+        real(kind(0d0)), intent(in) :: c
+        real(kind(0d0)) :: f_gauss_sigma_time_local
+        
+        if (gauss_conv_flag) then
+            f_gauss_sigma_time_local = gauss_sigma_dist(ai) / c
+        else
+            f_gauss_sigma_time_local = gauss_sigma_time(ai)
+        end if
+    end function f_gauss_sigma_time_local
+
 end module m_acoustic_src

toolchain/mfc/test/case.py

@@ -217,6 +217,9 @@ def compute_tolerance(self) -> float:
         if self.params.get("bubbles", 'F') == 'T':
             return 2e-10
 
+        if self.params.get("acoustic_source", 'F') == 'T':
+            return 2e-12
+
         if self.params.get("hypoelasticity", 'F') == 'T':
             return 1e-7