Surface Grid Refinement

examples/example_config_files/ivc_inactivelayer.yaml

@@ -20,6 +20,7 @@ grid:
       boundaries:
         left: inf
         right: inf
+  spacing_surface_refinement: [1e-4, 1e-4, 1e-4] #m
 medium: vacuum
 detectors:
 - semiconductor:

examples/example_config_files/public_CGD_config_cyl_grid.yaml

@@ -10,7 +10,7 @@ grid:
     r:
       to: 9
       boundaries: inf
-    y:
+    phi:
       from: 0
       to: 0
       boundaries: periodic

examples/example_config_files/test_detector.yaml

@@ -10,7 +10,7 @@ grid:
     r:
       to: 9
       boundaries: inf
-    y:
+    phi:
       from: 0
       to: 180
       boundaries: periodic

src/PotentialCalculation/Refinement.jl

@@ -174,3 +174,76 @@ end
 
 _extend_refinement_limits(rl::Real) = (rl, rl, rl )
 _extend_refinement_limits(rl::Tuple{<:Real,<:Real,<:Real}) = rl
+
+@inline function has_surface_points(slice::AbstractArray{PointType})::Bool
+    return any(is_in_inactive_layer, slice)
+end
+
+function _refine_axis_surface( ax::DiscreteAxis{T, <:Any, <:Any, ClosedInterval{T}}, surface_intervals::AbstractVector{Bool}, min_spacing::T;
+    extra_before::Int = 5,  # intervals to refine before first surface interval
+    extra_after::Int = 5    # intervals to refine after last surface interval
+) where {T}
+
+    old_ticks = ax.ticks
+    n_int = length(surface_intervals)
+
+    # Find first and last surface intervals
+    first_surface = findfirst(surface_intervals)
+    last_surface  = findlast(surface_intervals)
+
+    (isnothing(first_surface) || isnothing(last_surface)) && return ax
+
+    # Create flags for all intervals to refine
+    refine_flags = copy(surface_intervals)
+
+    # Add extra intervals before first surface interval
+    start_idx = max(first_surface - extra_before, 1)
+    refine_flags[start_idx:first_surface-1] .= true
+
+    # Add extra intervals after last surface interval
+    end_idx = min(last_surface + extra_after, n_int)
+    refine_flags[last_surface+1:end_idx] .= true
+
+    # Merge consecutive intervals to refine
+    merged = Vector{UnitRange{Int}}()
+    i = 1
+    while i <= n_int
+        if refine_flags[i]
+            j = i
+            while j < n_int && refine_flags[j+1]
+                j += 1
+            end
+            push!(merged, i:j)
+            i = j + 1
+        else
+            i += 1
+        end
+    end
+
+    # Compute number of points to add per interval
+    ns = zeros(Int, n_int)
+    for r in merged
+        for i in r
+            Δ = old_ticks[i+1] - old_ticks[i]
+            if Δ > min_spacing
+                ns[i] = ceil(Int, Δ/min_spacing) - 1
+            end
+        end
+    end
+
+    # Subdivide intervals
+    sub_widths = [(old_ticks[i+1] - old_ticks[i]) / (ns[i]+1) for i in 1:n_int]
+    ticks = Vector{T}(undef, length(old_ticks) + sum(ns))
+    i_tick = 1
+    for j in 1:n_int
+        ticks[i_tick] = old_ticks[j]
+        for k in 1:ns[j]
+            i_tick += 1
+            ticks[i_tick] = old_ticks[j] + k*sub_widths[j]
+        end
+        i_tick += 1
+    end
+    ticks[end] = old_ticks[end]
+
+    return typeof(ax)(ax.interval, ticks)
+end

src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl

@@ -10,14 +10,13 @@
     pwΔmp2r, pwΔmp2l, 
     pwΔmp3r, pwΔmp3l,
 ) where {T, S}
-
 
     # pww1r        = pcs.geom_weights[3][1, in1]
     pww1r        = geom_weights_3[1, in1]
     pww1l        = geom_weights_3[2, in1]
     pwΔmp1       = geom_weights_3[3, in1]
     Δ1_ext_inv_l = geom_weights_3[4, in1]
-    Δ1_ext_inv_r = geom_weights_3[4, in1 + 1]
+    Δ1_ext_inv_r = in1 < size(geom_weights_3, 2) ? geom_weights_3[4, in1 + 1] : geom_weights_3[4, in1]
 
     # ϵ_ijk: i: 1.RB-Dim. | j: 2.RB-Dim. | k: 3.RB-Dim.
     ϵ_lll, ϵ_llr, ϵ_lrl, ϵ_lrr, ϵ_rll, ϵ_rlr, ϵ_rrl, ϵ_rrr = get_ϵ_of_oktant(
@@ -168,36 +167,50 @@ end
     ϵ_r::AbstractArray{T, 3}, ::Type{Cylindrical},
     i1, in1, i2, in2, i3, in3
 ) where {T}
+    # Cylindrical: inner loop over z (3rd dimension)
+    i1_safe  = clamp(i1, 1, size(ϵ_r, 3))
+    in1_safe = clamp(in1, 1, size(ϵ_r, 3))
+    i2_safe  = clamp(i2, 1, size(ϵ_r, 2))
+    in2_safe = clamp(in2, 1, size(ϵ_r, 2))
+    i3_safe  = clamp(i3, 1, size(ϵ_r, 1))
+    in3_safe = clamp(in3, 1, size(ϵ_r, 1))
     # ϵ_r is not transformed into an red-black-4D-array.
     # The inner loop (over i1) is along the z-Dimension (Cylindrical Case), 
     # which is the 3rd dimension for Cylindrical coordinates: (r, φ, z)
     return @inbounds begin
-        ϵ_r[ in3, in2, in1 ],
-        ϵ_r[  i3, in2, in1 ],
-        ϵ_r[ in3,  i2, in1 ],
-        ϵ_r[  i3,  i2, in1 ],
-        ϵ_r[ in3, in2,  i1 ],
-        ϵ_r[  i3, in2,  i1 ],
-        ϵ_r[ in3,  i2,  i1 ],
-        ϵ_r[  i3,  i2,  i1 ]
+        ϵ_r[in3_safe, in2_safe, in1_safe],
+        ϵ_r[i3_safe,  in2_safe, in1_safe],
+        ϵ_r[in3_safe,  i2_safe, in1_safe],
+        ϵ_r[i3_safe,   i2_safe, in1_safe],
+        ϵ_r[in3_safe, in2_safe,  i1_safe],
+        ϵ_r[i3_safe,  in2_safe,  i1_safe],
+        ϵ_r[in3_safe,  i2_safe,  i1_safe],
+        ϵ_r[i3_safe,   i2_safe,  i1_safe]
     end
 end
 
 @inline function get_ϵ_of_oktant(
     ϵ_r::AbstractArray{T, 3}, ::Type{Cartesian},
     i1, in1, i2, in2, i3, in3
 ) where {T}
+    # Cartesian: inner loop over x (1st dimension)
+    i1_safe  = clamp(i1, 1, size(ϵ_r, 1))
+    in1_safe = clamp(in1, 1, size(ϵ_r, 1))
+    i2_safe  = clamp(i2, 1, size(ϵ_r, 2))
+    in2_safe = clamp(in2, 1, size(ϵ_r, 2))
+    i3_safe  = clamp(i3, 1, size(ϵ_r, 3))
+    in3_safe = clamp(in3, 1, size(ϵ_r, 3))
     # The inner loop (over i1) is along the x-Dimension (Cartesian Case), 
     # which is the 1rd dimension for Cartesian coordinates: (x, y, z)
     return @inbounds begin
-        ϵ_r[ in1, in2, in3 ], 
-        ϵ_r[ in1, in2,  i3 ],
-        ϵ_r[ in1,  i2, in3 ], 
-        ϵ_r[ in1,  i2,  i3 ],
-        ϵ_r[  i1, in2, in3 ],
-        ϵ_r[  i1, in2,  i3 ],
-        ϵ_r[  i1,  i2, in3 ],
-        ϵ_r[  i1,  i2,  i3 ]
+        ϵ_r[in1_safe, in2_safe, in3_safe],
+        ϵ_r[in1_safe, in2_safe,  i3_safe],
+        ϵ_r[in1_safe,  i2_safe, in3_safe],
+        ϵ_r[in1_safe,  i2_safe,  i3_safe],
+        ϵ_r[ i1_safe, in2_safe, in3_safe],
+        ϵ_r[ i1_safe, in2_safe,  i3_safe],
+        ϵ_r[ i1_safe,  i2_safe, in3_safe],
+        ϵ_r[ i1_safe,  i2_safe,  i3_safe]
     end
 end
 
@@ -242,4 +255,4 @@ end
     np::NTuple{6, T}, ::Type{Cartesian}, i::Int
 ) where {T}
     return np
-end
+end