A posible bug in xsmooth_sem?

[ZQiwen]

Hello everybody,

I am currently applying a Gaussian smoothing to the sensitivity kernels using the CPU version of xsmooth_sem compiled from codes in the directory src/tomography/postprocessing. I aim to suppress the singular values near the seismic stations (i.e. the adjoint sources) without too much blurring the kernels, so I choose small sigmas of 1000 meters in both v and h directions (which corresponds to a 2,800-meter smoothing radius) while the size of the element is ~5,000 meters. I may have found a bug that the GLL points from an adjacent element are probably not taken into account due to a misdefined element_size.

Here are the corresponding codes: Line 318-330 in smooth_sem.f90:

  ! adds margin to search radius
  element_size = max(sigma_h,sigma_v) * 0.5

  ! search radius
  sigma_h3 = 3.0  * sigma_h + element_size
  sigma_v3 = 3.0  * sigma_v + element_size

  ! helper variables
  sigma_h2_inv = 1.0_CUSTOM_REAL / sigma_h2
  sigma_v2_inv = 1.0_CUSTOM_REAL / sigma_v2

  sigma_h3_sq = sigma_h3 * sigma_h3
  sigma_v3_sq = sigma_v3 * sigma_v3

We can see that the element_size is defined to be half of the maximum input sigma rather than the size of an element, which seems to be confusing. The search radius defined here (3,500 meters in my settings) is then used to compare with the distance between the centers of two elements. Elements with distances larger than the search radius do not contribute to the smoothed kernel. See Line 898-903 in smooth_sem.f90:

! loop over slices do inum = 1,num_slices ... ! loop over elements to be smoothed in the current slice do ispec = 1, NSPEC_AB ... ! --- only double loop over the elements in the search radius --- do ispec2 = 1, NSPEC_N ...

          ! calculates horizontal and vertical distance between two element centers
          ! (squared distances)
          call get_distance_vec(dist_h,dist_v,center_x0,center_y0,center_z0,center_x,center_y,center_z)

          ! checks distance between centers of elements
          if (dist_h > sigma_h3_sq .or. dist_v > sigma_v3_sq) cycle

We can see that if two adjacent elements share a common face (the closest case), the distance between their centers is 5,000 meters in my case which is larger than the search radius (3,500 meters). Then the if-condition is satisfied and the element ispec2 is skipped from the later codes. This means the smoothed kernel receives no contribution even from the closest adjacent element although it should because the 2,800-meter smoothing radius is larger than the distance between two adjacent GLL points (~1,250 meters, NGLL=5).

This problem is rather only if the smoothing radius ( sqrt(8)*sigma ) is relatively small compared to the size of the element which is not common in real applications, or is it better to define the element_size to be the size of the element as its name is? For example:

  ! gets mesh dimensions
  call check_mesh_distances(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
                            x_min_glob,x_max_glob,y_min_glob,y_max_glob,z_min_glob,z_max_glob, &
                            elemsize_min_glob,elemsize_max_glob, &
                            distance_min_glob,distance_max_glob)
  call max_all_cr(elemsize_max_glob, element_size)

I hope I am making my point and do not hesitate to ask me if not. Cheers,

Qiwen

[danielpeter]

yes, that's a good point.

and correct, the current smoothing wouldn't look beyond the element if the input smoothing lengths are smaller than the element sizes. that would still leave a discontinuity in the smoothed output values when there is a velocity jump between 2 elements.

in your suggestion, you would want to smooth across elements on common GLL points. that makes sense. smoothing within single elements only is probably not helping much. let me see if we can just add the element size estimated by the check_mesh_distances routine to the search radius to solve this.

[ZQiwen]

Hi Daniel,

Thanks for responding.

I have made two mistakes in the previous comment, we should use

call max_all_all_cr(elemsize_max_glob*dsqrt(3.d0), element_size)

instead of

call max_all_cr(elemsize_max_glob, element_size)

use subroutine max_all_all_cr because the variable element_size is needed by all processors and times *dsqrt(3.d0) to elemsize_max_glob for cases when two elements share a common corner.

This is how it looks like when using a small smoothing length (10 meters, much smaller then the element size of 5000 meters) after modifying the definition of element_size as mentioned:

Eight values (black diamonds) at the same position (here by the Z coordinate) for the case when eight elements share a common corner and four values for sharing a common edge. We can see there are indeed jumps in values of kernel between two adjacent (share a common face/edge/corner) elements, which will not be smoothed out using the original smoothing code with such a small smoothing length.

[danielpeter]

thanks! your suggestion has been included in PR #1502

[ZQiwen]

Cheers!