diff --git a/docs/documentation/case.md b/docs/documentation/case.md index d7cdad48c2..5b166063b9 100644 --- a/docs/documentation/case.md +++ b/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. diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index e2fa6ff4e5..6d285bfc1e 100644 --- a/src/simulation/m_checker.fpp +++ b/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 '// & diff --git a/src/simulation/m_monopole.fpp b/src/simulation/m_monopole.fpp index 977da2c4b0..3acc14cd07 100644 --- a/src/simulation/m_monopole.fpp +++ b/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 diff --git a/toolchain/mfc/test/case.py b/toolchain/mfc/test/case.py index b49b15f3f0..7d2681f43b 100644 --- a/toolchain/mfc/test/case.py +++ b/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