Add gradient smoothing for elastic backend

src/models/elastic/backends/Elastic2D_Iso_CPML_AMDGPU.jl

@@ -11,5 +11,6 @@ using SeismicWaves.FiniteDifferencesMacros
 include("shared/standard_xPU.jl")
 include("shared/elastic2D_iso_xPU.jl")
 include("shared/correlate_gradient_xPU.jl")
+include("shared/smooth_gradient_2D.jl")
 
 end

src/models/elastic/backends/Elastic2D_Iso_CPML_CUDA.jl

@@ -11,5 +11,6 @@ using SeismicWaves.FiniteDifferencesMacros
 include("shared/standard_xPU.jl")
 include("shared/elastic2D_iso_xPU.jl")
 include("shared/correlate_gradient_xPU.jl")
+include("shared/smooth_gradient_2D.jl")
 
 end

src/models/elastic/backends/Elastic2D_Iso_CPML_Serial.jl

@@ -708,6 +708,8 @@ function correlate_gradients!(grid, vcurr, vold, dt, freetop)
     )
 end
 
+include("shared/smooth_gradient_2D.jl")
+
 #########################################
 end  # end module
 #########################################

src/models/elastic/backends/Elastic2D_Iso_CPML_Threads.jl

@@ -10,5 +10,6 @@ using SeismicWaves.FiniteDifferencesMacros
 include("shared/standard_xPU.jl")
 include("shared/elastic2D_iso_xPU.jl")
 include("shared/correlate_gradient_xPU.jl")
+include("shared/smooth_gradient_2D.jl")
 
 end

src/models/elastic/backends/shared/smooth_gradient_2D.jl

@@ -0,0 +1,22 @@
+function smooth_gradient!(grad, possrcs, radius)
+    nx, ny = size(grad)
+    nsources = size(possrcs, 1)
+    for s in 1:nsources
+        isrc = possrcs[s, 1]
+        jsrc = possrcs[s, 2]
+        imin = isrc - radius
+        imax = isrc + radius
+        jmin = jsrc - radius
+        jmax = jsrc + radius
+        for i in imin:imax
+            for j in jmin:jmax
+                if 1 <= i <= nx && 1 <= j <= ny
+                    r = sqrt((i - isrc)^2 + (j - jsrc)^2)
+                    if r < radius
+                        grad[i, j] *= r / radius
+                    end
+                end
+            end
+        end
+    end
+end

src/models/elastic/ela_gradient.jl

@@ -167,7 +167,11 @@ function swgradient_1shot!(
     gradient_ρ .+= back_interp(model.matprop.interp_method_ρ, model.matprop.ρ, Array(grid.fields["grad_ρ_ihalf"].value), 1) .+
                    back_interp(model.matprop.interp_method_ρ, model.matprop.ρ, Array(grid.fields["grad_ρ_jhalf"].value), 2)
     gradient_μ .+= back_interp(model.matprop.interp_method_μ, model.matprop.μ, Array(grid.fields["grad_μ_ihalf_jhalf"].value), [1, 2])
-    # TODO smooth gradients
+    # Smooth gradients
+    nearest_grdpts_src = find_nearest_grid_points(model, shot.srcs.positions)
+    backend.smooth_gradient!(gradient_ρ, nearest_grdpts_src, model.smooth_radius)
+    backend.smooth_gradient!(gradient_λ, nearest_grdpts_src, model.smooth_radius)
+    backend.smooth_gradient!(gradient_μ, nearest_grdpts_src, model.smooth_radius)
     # Compute regularization if needed
     ∂χ_∂ρ, ∂χ_∂λ, ∂χ_∂μ = (misfit.regularization !== nothing) ? dχ_dm(misfit.regularization, model.matprop) : (0, 0, 0)
     # Return gradients
@@ -336,7 +340,11 @@ function swgradient_1shot!(
     gradient_ρ .+= back_interp(model.matprop.interp_method_ρ, model.matprop.ρ, Array(grid.fields["grad_ρ_ihalf"].value), 1) .+
                    back_interp(model.matprop.interp_method_ρ, model.matprop.ρ, Array(grid.fields["grad_ρ_jhalf"].value), 2)
     gradient_μ .+= back_interp(model.matprop.interp_method_μ, model.matprop.μ, Array(grid.fields["grad_μ_ihalf_jhalf"].value), [1, 2])
-    # TODO smooth gradients
+    # Smooth gradients
+    nearest_grdpts_src = find_nearest_grid_points(model, shot.srcs.positions)
+    backend.smooth_gradient!(gradient_ρ, nearest_grdpts_src, model.smooth_radius)
+    backend.smooth_gradient!(gradient_λ, nearest_grdpts_src, model.smooth_radius)
+    backend.smooth_gradient!(gradient_μ, nearest_grdpts_src, model.smooth_radius)
     # Compute regularization if needed
     ∂χ_∂ρ, ∂χ_∂λ, ∂χ_∂μ = (misfit.regularization !== nothing) ? dχ_dm(misfit.regularization, model.matprop) : (0, 0, 0)
     # Return gradients