Add threshold and support width as inputs to the model patch and mono…

…pole parameters. (#336)

docs/documentation/case.md

@@ -175,6 +175,7 @@ To restart simulation from $k$-th time step, set `t_step_start`=k.
 | `model%rotate(i)`    | Real    | Not Supported         | Model's (applied) angle of rotation about axis $i$.          |
 | `model%translate(i)` | Real    | Not Supported         | Model's $i$-th component of (applied) translation.           |
 | `model%spc`          | Integer | Not Supported         | Number of samples per cell when discretizing the model into the grid. |
+| `model%threshold`    | Real    | Not Supported         | Ray fraction inside the model patch above which the fraction is set to one.|
 
 *: These parameters should be prepended with `patch_icpp(j)%` where $j$ is the patch index. 
 
@@ -451,17 +452,18 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 
 ### 7. Acoustic Source
 
-| Parameter         | Type    | Description |
-| ---:              | :----:  | :--- |
-| `Monopole`        | Logical | Acoustic source |
-| `num_mono`        | Integer | Number of acoustic sources |
-| `Mono(i)%pulse`   | Integer | Acoustic wave form: [1] Sine [2] Gaussian [3] Square |
-| `Mono(i)%npulse`  | Integer | Number of pulse cycles |
-| `Mono(i)%support` | Integer | Type of the spatial support of the acoustic source : [1] 1D [2] Finite width (2D) [3] Support for finite line/patch [4] General support for 3D simulation in cartesian systems [5] Support along monopole acoustic transducer [6] Support for cylindrical coordinate system along axial-dir |
-| `Mono(i)%loc(j)`  | Real    | $j$-th coordinate of the point that consists of $i$-th source plane |
-| `Mono(i)%dir`     | Real    | Direction of acoustic propagation	|
-| `Mono(i)%mag`     | Real    | Pulse magnitude	|
-| `Mono(i)%length`  | Real    | Spatial pulse length |
+| Parameter                | Type    | Description |
+| ---:                     | :----:  | :--- |
+| `Monopole`               | Logical | Acoustic source |
+| `num_mono`               | Integer | Number of acoustic sources |
+| `Mono(i)%pulse`          | Integer | Acoustic wave form: [1] Sine [2] Gaussian [3] Square |
+| `Mono(i)%npulse`         | Integer | Number of pulse cycles |
+| `Mono(i)%support`        | Integer | Type of the spatial support of the acoustic source : [1] 1D [2] Finite width (2D) [3] Support for finite line/patch [4] General support for 3D simulation in cartesian systems [5] Support along monopole acoustic transducer [6] Support for cylindrical coordinate system along axial-dir |
+| `Mono(i)%support_width`  | Real    | The width of the monopole support in terms of cell width |
+| `Mono(i)%loc(j)`         | Real    | $j$-th coordinate of the point that consists of $i$-th source plane |
+| `Mono(i)%dir`            | Real    | Direction of acoustic propagation	|
+| `Mono(i)%mag`            | Real    | Pulse magnitude	|
+| `Mono(i)%length`         | Real    | Spatial pulse length |
 
 The table lists acoustic source parameters.
 The parameters are optionally used to define a source plane in the domain that generates an acoustic wave that propagates in a specified direction normal to the source plane (one-way acoustic source).
@@ -490,6 +492,9 @@ The $i$-th source plane is determined by the point at [`Mono(i)%loc(1)`, `Mono(i
 The source plane is defined in the finite region of the domain: $x\in[-\infty,\infty]$ and $y,z\in$[-`mymono_length`/2, `mymono_length`/2].
 There are a few additional spatial support types available for special source types and coordinate systems tabulated in [Monopole supports](#monopole-supports).
 
+- `Mono(i)%support_width` defines how many cell width the monopole support function extended by.
+Large `Mono(i)%support_width` is preferred when `Mono(i)%mag` is large.
+
 ### 8. Ensemble-Averaged Bubble Model
 
 | Parameter      | Type    | Description                                    |

src/common/m_derived_types.fpp

@@ -94,6 +94,9 @@ module m_derived_types
 
         integer :: spc !<
         !! Number of samples per cell to use when discretizing the STL object.
+
+        real(kind(0d0)) :: threshold !<
+        !! Threshold to turn on smoothen STL patch.
     end type ic_model_parameters
 
     type :: t_triangle
@@ -261,6 +264,7 @@ module m_derived_types
         integer :: support
         real(kind(0d0)) :: aperture
         real(kind(0d0)) :: foc_length
+        real(kind(0d0)) :: support_width
     end type mono_parameters
 
     !> Ghost Point for Immersed Boundaries

src/pre_process/m_global_parameters.fpp

@@ -303,6 +303,7 @@ contains
             patch_icpp(i)%model%translate(:) = 0d0
             patch_icpp(i)%model%filepath(:) = ' '
             patch_icpp(i)%model%spc = 10
+            patch_icpp(i)%model%threshold = 0.9d0
             patch_icpp(i)%x_centroid = dflt_real
             patch_icpp(i)%y_centroid = dflt_real
             patch_icpp(i)%z_centroid = dflt_real

src/pre_process/m_mpi_proxy.fpp

@@ -65,7 +65,7 @@ contains
             & 'bc_x%end', 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end',      &
             & 'perturb_flow_mag', 'pref', 'rhoref', 'poly_sigma', 'R0ref',     &
             & 'Web', 'Ca', 'Re_inv', 'sigR', 'sigV', 'rhoRV', 'palpha_eps',    &
-            & 'ptgalpha_eps' ]
+            & 'ptgalpha_eps']
             call MPI_BCAST(${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
@@ -80,7 +80,8 @@ contains
             #:for VAR in [ 'x_centroid', 'y_centroid', 'z_centroid',           &
                 & 'length_x', 'length_y', 'length_z', 'radius', 'epsilon',     &
                 & 'beta', 'smooth_coeff', 'rho', 'p0', 'm0', 'r0', 'v0',       &
-                & 'pres', 'gamma', 'pi_inf', 'hcid', 'cv', 'qv', 'qvp' ]
+                & 'pres', 'gamma', 'pi_inf', 'hcid', 'cv', 'qv', 'qvp',        &
+                & 'model%threshold']
                 call MPI_BCAST(patch_icpp(i)%${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 

src/pre_process/m_patches.fpp

@@ -1942,6 +1942,18 @@ contains
 
                     eta = f_model_is_inside(model, point, (/dx, dy, dz/), patch_icpp(patch_id)%model%spc)
 
+                    if (patch_icpp(patch_id)%smoothen) then
+                        if (eta > patch_icpp(patch_id)%model%threshold) then
+                            eta = 1d0
+                        end if
+                    else
+                        if (eta > patch_icpp(patch_id)%model%threshold) then
+                            eta = 1d0
+                        else
+                            eta = 0d0
+                        end if
+                    end if
+
                     call s_assign_patch_primitive_variables(patch_id, i, j, k, &
                                                             eta, q_prim_vf, patch_id_fp)
 

src/simulation/m_global_parameters.fpp

@@ -491,6 +491,10 @@ contains
             mono(j)%support = 1
             mono(j)%foc_length = dflt_real
             mono(j)%aperture = dflt_real
+            ! The author suggested the support width is typically on the order of
+            ! the width of the characteristic cells. Here, we choose 2.5 cell width
+            ! as the default value.
+            mono(j)%support_width = 2.5d0
         end do
 
         fd_order = dflt_int

src/simulation/m_monopole.fpp

@@ -25,8 +25,8 @@ module m_monopole
     real(kind(0d0)), allocatable, dimension(:, :) :: loc_mono
     !$acc declare create(loc_mono)
 
-    real(kind(0d0)), allocatable, dimension(:) :: foc_length, aperture
-    !$acc declare create(foc_length, aperture)
+    real(kind(0d0)), allocatable, dimension(:) :: foc_length, aperture, support_width
+    !$acc declare create(foc_length, aperture, support_width)
 
     real(kind(0d0)), allocatable, dimension(:) :: mag, length, npulse, dir, delay
     !$acc declare create(mag, length, npulse, dir, delay)
@@ -36,7 +36,7 @@ contains
     subroutine s_initialize_monopole_module()
         integer :: i, j !< generic loop variables
 
-        @:ALLOCATE(mag(1:num_mono), support(1:num_mono), length(1:num_mono), npulse(1:num_mono), pulse(1:num_mono), dir(1:num_mono), delay(1:num_mono), loc_mono(1:3, 1:num_mono), foc_length(1:num_mono), aperture(1:num_mono))
+        @:ALLOCATE(mag(1:num_mono), support(1:num_mono), length(1:num_mono), npulse(1:num_mono), pulse(1:num_mono), dir(1:num_mono), delay(1:num_mono), loc_mono(1:3, 1:num_mono), foc_length(1:num_mono), aperture(1:num_mono), support_width(1:num_mono))
 
         do i = 1, num_mono
             mag(i) = mono(i)%mag
@@ -48,11 +48,12 @@ contains
             delay(i) = mono(i)%delay
             foc_length(i) = mono(i)%foc_length
             aperture(i) = mono(i)%aperture
+            support_width(i) = mono(i)%support_width
             do j = 1, 3
                 loc_mono(j, i) = mono(i)%loc(j)
             end do
         end do
-        !$acc update device(mag, support, length, npulse, pulse, dir, delay, foc_length, aperture, loc_mono)
+        !$acc update device(mag, support, length, npulse, pulse, dir, delay, foc_length, aperture, loc_mono, support_width)
 
     end subroutine
 
@@ -327,13 +328,13 @@ contains
 
         if (n == 0) then
             sig = dx(j)
-            sig = sig*2.5d0
+            sig = sig*support_width(nm)
         else if (p == 0) then
             sig = maxval((/dx(j), dy(k)/))
-            sig = sig*2.5d0
+            sig = sig*support_width(nm)
         else
             sig = maxval((/dx(j), dy(k), dz(l)/))
-            sig = sig*2.5d0
+            sig = sig*support_width(nm)
         end if
 
         if (n == 0) then      !1D
@@ -441,7 +442,7 @@ contains
                 ! Support for cylindrical coordinate system
                 !sig = maxval((/dx(j), dy(k)*sin(dz(l)), dz(l)*cos(dz(l))/))
                 sig = dx(j)
-                sig = sig*2.5d0
+                sig = sig*support_width(nm)
                 hx_cyl = x_cc(j) - mono_loc(1)
                 hy_cyl = y_cc(k)*sin(z_cc(l)) - mono_loc(2)
                 hz_cyl = y_cc(k)*cos(z_cc(l)) - mono_loc(3)

src/simulation/m_mpi_proxy.fpp

@@ -184,7 +184,7 @@ contains
             end do
 
             #:for VAR in [ 'mag', 'length', 'delay', 'dir', 'npulse', 'pulse',  &
-                'support', 'foc_length', 'aperture' ]
+                'support', 'foc_length', 'aperture', 'support_width' ]
                 call MPI_BCAST(mono(j)%${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 

toolchain/mfc/run/case_dicts.py

@@ -58,6 +58,7 @@
             PRE_PROCESS.append(f"patch_icpp({p_id})%model%{attribute}({j})")
 
     PRE_PROCESS.append(f"patch_icpp({p_id})%model%spc")
+    PRE_PROCESS.append(f"patch_icpp({p_id})%model%threshold")
 
     for cmp_id, cmp in enumerate(["x", "y", "z"]):
         cmp_id += 1
@@ -128,7 +129,7 @@
 
     for mono_id in range(1,4+1):
         for attribute in ["mag", "length", "dir", "npulse", "pulse", "support",
-                          "delay", "foc_length", "aperture"]:
+                          "delay", "foc_length", "aperture", "support_width"]:
             SIMULATION.append(f"Mono({mono_id})%{attribute}")
 
         for cmp_id in range(1,3+1):