From d864051d4f9ff210c894d063004b53f4c0584b78 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 17 Feb 2026 13:00:40 +0100
Subject: [PATCH 01/12] Control surface refinement in yaml

---
 examples/example_config_files/ivc_inactivelayer.yaml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/example_config_files/ivc_inactivelayer.yaml b/examples/example_config_files/ivc_inactivelayer.yaml
index b99ee6378..9336aecfe 100644
--- a/examples/example_config_files/ivc_inactivelayer.yaml
+++ b/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]
 medium: vacuum
 detectors:
 - semiconductor:

From 44e7dc1fe14cf5637edb13e76e52de02ec5e6038 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 17 Feb 2026 13:01:13 +0100
Subject: [PATCH 02/12] Detect and refine surface intervals

---
 src/PotentialCalculation/Refinement.jl | 71 ++++++++++++++++++++++++++
 1 file changed, 71 insertions(+)

diff --git a/src/PotentialCalculation/Refinement.jl b/src/PotentialCalculation/Refinement.jl
index fab5dac6b..e8ca9f8e0 100644
--- a/src/PotentialCalculation/Refinement.jl
+++ b/src/PotentialCalculation/Refinement.jl
@@ -174,3 +174,74 @@ 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) where {T}
+
+    old_ticks = ax.ticks
+    n_int = length(surface_intervals)
+
+    # Merge consecutive surface intervals
+    merged = Vector{UnitRange{Int}}()
+    i = 1
+    while i <= n_int
+        if surface_intervals[i]
+            j = i
+            while j < n_int && surface_intervals[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
+
+    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
+
+function _create_refined_grid_surface(p::ElectricPotential, point_types::Array{PointType,3}, min_spacing::NTuple{3,T}) where T
+
+    sz = size(point_types)
+    surface_intervals = ntuple(d -> falses(sz[d]-1), 3)
+
+    for i in 1:sz[1]-1
+        surface_intervals[1][i] = has_surface_points(point_types[i:i+1, :, :])
+    end
+    for i in 1:sz[2]-1
+        surface_intervals[2][i] = has_surface_points(point_types[:, i:i+1, :])
+    end
+    for i in 1:sz[3]-1
+        surface_intervals[3][i] = has_surface_points(point_types[:, :, i:i+1])
+    end
+    new_axes = ntuple(i -> _refine_axis_surface(p.grid.axes[i], surface_intervals[i], min_spacing[i]), 3)
+
+    return typeof(p.grid)(new_axes)
+end

From c91b326288e00b868204d82e556a08e309d5a9d7 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 17 Feb 2026 13:01:49 +0100
Subject: [PATCH 03/12] Take into account boundaries for small refinements
 1e-5m

---
 .../SuccessiveOverRelaxation/CPU_innerloop.jl | 51 ++++++++++++-------
 1 file changed, 32 insertions(+), 19 deletions(-)

diff --git a/src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl b/src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl
index e6bc1476c..fcb906309 100644
--- a/src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl
+++ b/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,18 +167,25 @@ 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
 
@@ -187,17 +193,24 @@ end
     ϵ_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
\ No newline at end of file
+end

From 15525311a78a418936e457403d573b8499cef6c9 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 17 Feb 2026 13:02:08 +0100
Subject: [PATCH 04/12] Handle surface refinement

---
 src/Simulation/Simulation.jl | 126 ++++++++++++++++++++++++++++++++++-
 1 file changed, 123 insertions(+), 3 deletions(-)

diff --git a/src/Simulation/Simulation.jl b/src/Simulation/Simulation.jl
index a8cf25af9..f3b75b72f 100644
--- a/src/Simulation/Simulation.jl
+++ b/src/Simulation/Simulation.jl
@@ -38,6 +38,7 @@ mutable struct Simulation{T <: SSDFloat, CS <: AbstractCoordinateSystem} <: Abst
     electric_potential::Union{ElectricPotential{T}, Missing}
     weighting_potentials::Vector{Any}
     electric_field::Union{ElectricField{T}, Missing}
+    spacing_surface_refinement::Union{NTuple{3,T}, Nothing}
 end
 
 function Simulation{T,CS}() where {T <: SSDFloat, CS <: AbstractCoordinateSystem}
@@ -54,7 +55,8 @@ function Simulation{T,CS}() where {T <: SSDFloat, CS <: AbstractCoordinateSystem
         missing,
         missing,
         [missing],
-        missing
+        missing,
+        nothing
     )
 end
 
@@ -113,6 +115,7 @@ function Simulation(nt::NamedTuple)
     else
         [missing for contact in sim.detector.contacts]
     end
+    sim.spacing_surface_refinement = get(nt, :spacing_surface_refinement, nothing)
     return sim
 end
 Base.convert(T::Type{Simulation}, x::NamedTuple) = T(x)
@@ -178,13 +181,19 @@ function Simulation{T}(dict::AbstractDict)::Simulation{T} where {T <: SSDFloat}
     sim.medium = material_properties[materials[haskey(dict, "medium") ? dict["medium"] : "vacuum"]]
     sim.detector = SolidStateDetector{T}(dict, sim.input_units) 
     sim.world = if haskey(dict, "grid") && isa(dict["grid"], AbstractDict) && haskey(dict["grid"], "axes")
-            World(T, dict["grid"], sim.input_units)
+        World(T, dict["grid"], sim.input_units)
         else let det = sim.detector 
             world_limits = get_world_limits_from_objects(CS, det)
             World(CS, world_limits)
         end
     end
     sim.weighting_potentials = Missing[ missing for i in 1:length(sim.detector.contacts)]
+    sim.spacing_surface_refinement = if haskey(dict, "grid") && isa(dict["grid"], AbstractDict) &&
+                                       haskey(dict["grid"], "spacing_surface_refinement")
+        ntuple(i -> T(dict["grid"]["spacing_surface_refinement"][i]), 3)
+    else
+        nothing
+    end
     return sim
 end
 
@@ -759,6 +768,7 @@ Takes the current state of `sim.electric_potential` and refines it with respect
 SolidStateDetectors.refine!(sim, ElectricPotential, max_diffs = (100, 100, 100), minimum_distances = (0.01, 0.02, 0.01))
 ```
 """
+
 function refine!(sim::Simulation{T}, ::Type{ElectricPotential},
                     max_diffs::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0)),
                     minimum_distances::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0));
@@ -771,6 +781,7 @@ function refine!(sim::Simulation{T}, ::Type{ElectricPotential},
         pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data,
                                         not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts)
 
+        
         sim.imp_scale = ImpurityScale(ImpurityScaleArray(pcs), sim.electric_potential.grid)
         sim.q_eff_imp = EffectiveChargeDensity(EffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
         sim.q_eff_fix = EffectiveChargeDensity(FixedEffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
@@ -812,6 +823,82 @@ function refine!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int
     nothing
 end
 
+"""
+Refine the simulation potential on the surface only
+"""
+
+function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4));
+                    not_only_paint_contacts::Bool = true,
+                    paint_contacts::Bool = true,
+                    update_other_fields::Bool = true) where {T <: SSDFloat}
+
+    old_grid = sim.electric_potential.grid
+    #Enforce lower bound of 1e-5
+    for d in 1:3
+        if length(old_grid.axes[d].ticks) != 1 && min_spacing[d] < T(1e-5)
+            @warn "min_spacing[$d] = $(min_spacing[d]) < 1e-5; uaing 1e-5m spacing."
+        end
+    end
+    min_spacing = ntuple(d -> length(old_grid.axes[d].ticks) == 1 ? T(0.0) : max(T(min_spacing[d]), T(1e-5)), 3)
+    
+    # Determine which intervals have surface points
+    surface_intervals = ntuple(d -> begin
+        n_int = length(old_grid.axes[d].ticks) - 1
+        map(i -> begin
+            slice = if d == 1
+                view(sim.point_types.data, i, :, :)
+            elseif d == 2
+                view(sim.point_types.data, :, i, :)
+            else
+                view(sim.point_types.data, :, :, i)
+            end
+            has_surface_points(slice)
+        end, 1:n_int)
+    end, 3)
+
+    # Refine axes
+    new_axes = ntuple(i -> _refine_axis_surface(old_grid.axes[i], surface_intervals[i], min_spacing[i]), 3)
+    new_grid = typeof(old_grid)(new_axes)
+
+    # Interpolate potential onto new grid
+    closed_pot = _get_closed_potential(sim.electric_potential)
+    new_data = Array{T,3}(undef, size(new_grid))
+    only2d = size(closed_pot.data, 2) == 1
+    int = interpolate_closed_potential(closed_pot, Val(only2d))
+
+    if only2d
+        for i3 in axes(new_data,3), i1 in axes(new_data,1)
+            new_data[i1,1,i3] = int(new_grid.axes[1].ticks[i1], new_grid.axes[3].ticks[i3])
+        end
+    else
+        for i3 in axes(new_data,3), i2 in axes(new_data,2), i1 in axes(new_data,1)
+            new_data[i1,i2,i3] = int(new_grid.axes[1].ticks[i1],
+                                      new_grid.axes[2].ticks[i2],
+                                      new_grid.axes[3].ticks[i3])
+        end
+    end
+
+    # Update electric potential in simulation
+    sim.electric_potential = _convert_to_original_potential(sim.electric_potential, new_data, new_grid)
+    
+    # Update dependent fields
+    if update_other_fields
+        
+        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data,
+                                        not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts)
+        sim.imp_scale = ImpurityScale(ImpurityScaleArray(pcs), sim.electric_potential.grid)
+        sim.q_eff_imp = EffectiveChargeDensity(EffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
+        sim.q_eff_fix = EffectiveChargeDensity(FixedEffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
+        sim.ϵ_r = DielectricDistribution(DielectricDistributionArray(pcs),
+                                        get_extended_midpoints_grid(sim.electric_potential.grid))
+        sim.point_types = PointTypes(PointTypeArray(pcs), sim.electric_potential.grid)
+    end
+
+    return nothing
+end
+
+
+
 """
     compute_min_tick_distance(grid)
 
@@ -899,6 +986,15 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
         else
             sor_consts = T.(sor_consts)
         end
+        has_surface_model = isdefined(sim.detector.semiconductor.impurity_density_model, :surface_imp_model)
+        if has_surface_model && (length(refinement_limits) < 3 || last(refinement_limits) > 0.05)
+            @warn """Surface model detected:
+                - Number of refinement steps: $(length(refinement_limits))
+                - Last refinement: $(last(refinement_limits))
+                        Surface refinement requires at least 3 refinement steps and the last refinement ≤ 0.05."""
+            refinement_limits = [0.2, 0.1, 0.05]
+            @warn """Falling back to default refinement_limits = $(refinement_limits)"""
+        end
       
         new_min_tick_distance::NTuple{3,T} = begin
             if ismissing(min_tick_distance)
@@ -978,6 +1074,9 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                                     device_array_type = device_array_type,
                                     use_nthreads = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.electric_potential.grid), max_nthreads[1], CS) : max_nthreads[1],
                                     sor_consts = sor_consts )
+
+
+            
         else
             apply_initial_state!(sim, potential_type, contact_id, grid; not_only_paint_contacts, paint_contacts, depletion_handling)
             update_till_convergence!( sim, potential_type, contact_id, convergence_limit,
@@ -1003,6 +1102,7 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                 else
                     abs.(ref_limits .* bias_voltage)
                 end
+                
                 refine!(sim, ElectricPotential, max_diffs, new_min_tick_distance)
                 nt = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.electric_potential.grid), max_nthreads[iref+1], CS) : max_nthreads[iref+1]
                 verbose && println("Grid size: $(size(sim.electric_potential.data)) - $(onCPU ? "using $(nt) threads now" : "GPU")") 
@@ -1015,6 +1115,26 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                                                 not_only_paint_contacts = not_only_paint_contacts, 
                                                 paint_contacts = paint_contacts,
                                                 sor_consts = is_last_ref ? T(1) : sor_consts )
+                
+                if has_surface_model && iref==3 
+                    mark_bulk_bits!(sim.point_types.data)
+                    mark_undep_bits!(sim.point_types.data, sim.imp_scale.data)
+                    mark_inactivelayer_bits!(sim.point_types.data)
+                    verbose && println("Surface Refinement")
+                    # Maximum spacing between (refined) surface ticks (if min_spacing = 1e-4m, only surface intervals wider than 0.1mm get new ticks)
+                    if sim.spacing_surface_refinement !== nothing
+                        refine_surface!(sim, sim.spacing_surface_refinement; update_other_fields=true)
+                    else
+                        refine_surface!(sim; update_other_fields=true)
+                    end
+                    update_till_convergence!( sim, potential_type, convergence_limit,
+                                              n_iterations_between_checks = n_iterations_between_checks,
+                                              max_n_iterations = max_n_iterations,
+                                              depletion_handling = depletion_handling,
+                                              device_array_type = device_array_type,
+                                              use_nthreads = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.electric_potential.grid),
+                                              max_nthreads[1], CS) : max_nthreads[1], sor_consts = T(1) )
+                end
             else
                 max_diffs = abs.(ref_limits)
                 refine!(sim, WeightingPotential, contact_id, max_diffs, new_min_tick_distance)
@@ -1057,7 +1177,7 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
     if depletion_handling && isEP
         mark_undep_bits!(sim.point_types.data, sim.imp_scale.data)
         
-        if isdefined(sim.detector.semiconductor.impurity_density_model, :surface_imp_model)
+        if has_surface_model
             mark_inactivelayer_bits!(sim.point_types.data)
         end
     end

From 4eae2a259126091ced09586ae5cd24693852338a Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 19 Feb 2026 14:22:28 +0100
Subject: [PATCH 05/12] Clean code and surface surrounding refinement

---
 src/PotentialCalculation/Refinement.jl | 56 ++++++++++++++------------
 1 file changed, 30 insertions(+), 26 deletions(-)

diff --git a/src/PotentialCalculation/Refinement.jl b/src/PotentialCalculation/Refinement.jl
index e8ca9f8e0..055ed0a7f 100644
--- a/src/PotentialCalculation/Refinement.jl
+++ b/src/PotentialCalculation/Refinement.jl
@@ -179,22 +179,44 @@ _extend_refinement_limits(rl::Tuple{<:Real,<:Real,<:Real}) = rl
     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) where {T}
+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)
 
-    # Merge consecutive surface intervals
+    # Find first and last surface intervals
+    first_surface = findfirst(surface_intervals)
+    last_surface  = findlast(surface_intervals)
+
+    if first_surface === nothing || last_surface === nothing
+        return ax # Nothing to refine
+    end
+
+    # 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 surface_intervals[i]
+        if refine_flags[i]
             j = i
-            while j < n_int && surface_intervals[j+1]
+            while j < n_int && refine_flags[j+1]
                 j += 1
             end
             push!(merged, i:j)
-            i = j+1
+            i = j + 1
         else
             i += 1
         end
@@ -211,8 +233,8 @@ function _refine_axis_surface(ax::DiscreteAxis{T, <:Any, <:Any, ClosedInterval{T
         end
     end
 
-    sub_widths = [(old_ticks[i+1]-old_ticks[i]) / (ns[i]+1) for i in 1:n_int]
-
+    # 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
@@ -224,24 +246,6 @@ function _refine_axis_surface(ax::DiscreteAxis{T, <:Any, <:Any, ClosedInterval{T
         i_tick += 1
     end
     ticks[end] = old_ticks[end]
-    return typeof(ax)(ax.interval, ticks)
-end
 
-function _create_refined_grid_surface(p::ElectricPotential, point_types::Array{PointType,3}, min_spacing::NTuple{3,T}) where T
-
-    sz = size(point_types)
-    surface_intervals = ntuple(d -> falses(sz[d]-1), 3)
-
-    for i in 1:sz[1]-1
-        surface_intervals[1][i] = has_surface_points(point_types[i:i+1, :, :])
-    end
-    for i in 1:sz[2]-1
-        surface_intervals[2][i] = has_surface_points(point_types[:, i:i+1, :])
-    end
-    for i in 1:sz[3]-1
-        surface_intervals[3][i] = has_surface_points(point_types[:, :, i:i+1])
-    end
-    new_axes = ntuple(i -> _refine_axis_surface(p.grid.axes[i], surface_intervals[i], min_spacing[i]), 3)
-
-    return typeof(p.grid)(new_axes)
+    return typeof(ax)(ax.interval, ticks)
 end

From 02b46dd8622744c83da616ac26385845ff06323c Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 24 Feb 2026 15:03:58 +0100
Subject: [PATCH 06/12] Fixed yaml files

---
 examples/example_config_files/ivc_inactivelayer.yaml          | 2 +-
 examples/example_config_files/public_CGD_config_cyl_grid.yaml | 2 +-
 examples/example_config_files/test_detector.yaml              | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/examples/example_config_files/ivc_inactivelayer.yaml b/examples/example_config_files/ivc_inactivelayer.yaml
index 9336aecfe..3bc9aff87 100644
--- a/examples/example_config_files/ivc_inactivelayer.yaml
+++ b/examples/example_config_files/ivc_inactivelayer.yaml
@@ -20,7 +20,7 @@ grid:
       boundaries:
         left: inf
         right: inf
-  spacing_surface_refinement: [1e-4, 1e-4, 1e-4]
+  spacing_surface_refinement: [1e-4, 1e-4, 1e-4] #m
 medium: vacuum
 detectors:
 - semiconductor:
diff --git a/examples/example_config_files/public_CGD_config_cyl_grid.yaml b/examples/example_config_files/public_CGD_config_cyl_grid.yaml
index 999663a1b..c03767cac 100644
--- a/examples/example_config_files/public_CGD_config_cyl_grid.yaml
+++ b/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
diff --git a/examples/example_config_files/test_detector.yaml b/examples/example_config_files/test_detector.yaml
index b50c18d08..67d2783e2 100644
--- a/examples/example_config_files/test_detector.yaml
+++ b/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

From cb008178a32232d7e897b6ea6a72e8098a5dc5de Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 24 Feb 2026 15:04:25 +0100
Subject: [PATCH 07/12] Support partial phi defined detectors

---
 src/Simulation/Simulation.jl | 116 +++++++++++++++++++++++++----------
 1 file changed, 84 insertions(+), 32 deletions(-)

diff --git a/src/Simulation/Simulation.jl b/src/Simulation/Simulation.jl
index f3b75b72f..000a5c6bf 100644
--- a/src/Simulation/Simulation.jl
+++ b/src/Simulation/Simulation.jl
@@ -824,23 +824,54 @@ function refine!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int
 end
 
 """
-Refine the simulation potential on the surface only
-"""
-
-function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4));
+refine_surface!(sim::Simulation{T};
+                    min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4)),
                     not_only_paint_contacts::Bool = true,
                     paint_contacts::Bool = true,
                     update_other_fields::Bool = true) where {T <: SSDFloat}
 
+Refine the simulation's electric potential grid **only in regions containing surface points**, keeping other regions unchanged.
+
+This function:
+
+1. Identifies which intervals along each axis contain surface points.
+2. Refines these intervals by adding extra grid points while respecting `min_spacing`.
+3. Interpolates the existing electric potential onto the refined grid.
+4. Updates the simulation's `electric_potential` and optionally dependent fields (`imp_scale`, `q_eff_imp`, `q_eff_fix`, `ϵ_r`, `point_types`).
+
+Special behaviour:
+
+- In cylindrical coordinates, the φ-axis (2nd axis) is left unchanged except for closing the axis for interpolation purposes.
+- Handles both 2D and 3D grids.
+
+# Arguments
+- `sim::Simulation{T}`: The simulation object containing the electric potential to refine.
+- `min_spacing::NTuple{3,T}`: Minimum allowed spacing in the surface after refinement. Default is `(1e-4, 1e-4, 1e-4)m`.
+"""
+
+function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4));
+                         not_only_paint_contacts::Bool = true,
+                         paint_contacts::Bool = true,
+                         update_other_fields::Bool = true) where {T <: SSDFloat}
+    
     old_grid = sim.electric_potential.grid
-    #Enforce lower bound of 1e-5
+
+    # Enforce minimum and maximum spacing for surface refinement (10 µm < min_spacing < 1 mm)
     for d in 1:3
-        if length(old_grid.axes[d].ticks) != 1 && min_spacing[d] < T(1e-5)
-            @warn "min_spacing[$d] = $(min_spacing[d]) < 1e-5; uaing 1e-5m spacing."
+        if length(old_grid.axes[d].ticks) != 1
+            if min_spacing[d] < T(1e-5) || min_spacing[d] > T(1e-3)
+                @warn "min_spacing[$d] = $(min_spacing[d]) m is out of bounds (1e-5 m < min_spacing < 1e-3 m). Higher values don't require surface refinement"
+                if min_spacing[d] < T(1e-5)
+                    @warn "Using min_spacing[$d] = 1e-5 m"
+                else
+                    @warn "Using min_spacing[$d] = 1e-3 m"  
+                end
+            end
         end
     end
-    min_spacing = ntuple(d -> length(old_grid.axes[d].ticks) == 1 ? T(0.0) : max(T(min_spacing[d]), T(1e-5)), 3)
     
+    min_spacing = ntuple(d -> length(old_grid.axes[d].ticks) == 1 ? T(0.0) : clamp(T(min_spacing[d]), T(1e-5), T(1e-3)), 3)
+
     # Determine which intervals have surface points
     surface_intervals = ntuple(d -> begin
         n_int = length(old_grid.axes[d].ticks) - 1
@@ -856,41 +887,55 @@ function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4)
         end, 1:n_int)
     end, 3)
 
-    # Refine axes
-    new_axes = ntuple(i -> _refine_axis_surface(old_grid.axes[i], surface_intervals[i], min_spacing[i]), 3)
-    new_grid = typeof(old_grid)(new_axes)
+    # Refine axes: skip phi-axis if Cylindrical
+    is_cyl = old_grid isa Grid{<:Any, 3, Cylindrical}
+    new_axes = ntuple(i -> begin
+                          if is_cyl && i == 2
+                              # Keep φ-axis type but close it for interpolation
+                              _get_closed_axis(old_grid.axes[i])
+                          else
+                              _refine_axis_surface( old_grid.axes[i], surface_intervals[i], min_spacing[i],)
+                          end
+                      end, 3)
+    
+    if is_cyl
+        new_grid = Grid{T,3,Cylindrical}((new_axes[1], new_axes[2], new_axes[3]))
+    else
+        new_grid = Grid{T,3,Cartesian}((new_axes[1], new_axes[2], new_axes[3]))
+    end
 
     # Interpolate potential onto new grid
     closed_pot = _get_closed_potential(sim.electric_potential)
     new_data = Array{T,3}(undef, size(new_grid))
     only2d = size(closed_pot.data, 2) == 1
     int = interpolate_closed_potential(closed_pot, Val(only2d))
-
+    
     if only2d
+        # 2D grid
         for i3 in axes(new_data,3), i1 in axes(new_data,1)
             new_data[i1,1,i3] = int(new_grid.axes[1].ticks[i1], new_grid.axes[3].ticks[i3])
         end
     else
+        # 3D grid
         for i3 in axes(new_data,3), i2 in axes(new_data,2), i1 in axes(new_data,1)
             new_data[i1,i2,i3] = int(new_grid.axes[1].ticks[i1],
-                                      new_grid.axes[2].ticks[i2],
-                                      new_grid.axes[3].ticks[i3])
+                                     new_grid.axes[2].ticks[i2],
+                                     new_grid.axes[3].ticks[i3])
         end
     end
-
+    
     # Update electric potential in simulation
     sim.electric_potential = _convert_to_original_potential(sim.electric_potential, new_data, new_grid)
     
     # Update dependent fields
     if update_other_fields
-        
+
         pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data,
                                         not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts)
         sim.imp_scale = ImpurityScale(ImpurityScaleArray(pcs), sim.electric_potential.grid)
         sim.q_eff_imp = EffectiveChargeDensity(EffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
         sim.q_eff_fix = EffectiveChargeDensity(FixedEffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
-        sim.ϵ_r = DielectricDistribution(DielectricDistributionArray(pcs),
-                                        get_extended_midpoints_grid(sim.electric_potential.grid))
+        sim.ϵ_r = DielectricDistribution(DielectricDistributionArray(pcs), get_extended_midpoints_grid(sim.electric_potential.grid))
         sim.point_types = PointTypes(PointTypeArray(pcs), sim.electric_potential.grid)
     end
 
@@ -987,15 +1032,25 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
             sor_consts = T.(sor_consts)
         end
         has_surface_model = isdefined(sim.detector.semiconductor.impurity_density_model, :surface_imp_model)
-        if has_surface_model && (length(refinement_limits) < 3 || last(refinement_limits) > 0.05)
-            @warn """Surface model detected:
-                - Number of refinement steps: $(length(refinement_limits))
-                - Last refinement: $(last(refinement_limits))
-                        Surface refinement requires at least 3 refinement steps and the last refinement ≤ 0.05."""
-            refinement_limits = [0.2, 0.1, 0.05]
-            @warn """Falling back to default refinement_limits = $(refinement_limits)"""
+
+        if has_surface_model    
+            if sim.spacing_surface_refinement === nothing
+                @warn """Surface model detected but `spacing_surface_refinement` is not defined.
+                         Surface refinement will not be performed."""
+            else
+                if length(refinement_limits) < 3 || last(refinement_limits) > 0.05
+                    @warn """Surface model detected:
+                             - Number of refinement steps: $(length(refinement_limits))
+                             - Last refinement: $(last(refinement_limits))
+                             Surface refinement requires:
+                             - At least 3 refinement steps
+                             - Last refinement ≤ 0.05"""
+                    refinement_limits = [0.2, 0.1, 0.05]
+                    @warn "Falling back to default refinement_limits = $(refinement_limits)"
+                end
+            end    
         end
-      
+        
         new_min_tick_distance::NTuple{3,T} = begin
             if ismissing(min_tick_distance)
                 compute_min_tick_distance(grid)
@@ -1116,17 +1171,14 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                                                 paint_contacts = paint_contacts,
                                                 sor_consts = is_last_ref ? T(1) : sor_consts )
                 
-                if has_surface_model && iref==3 
+                if has_surface_model && iref==3 && sim.spacing_surface_refinement !== nothing
                     mark_bulk_bits!(sim.point_types.data)
                     mark_undep_bits!(sim.point_types.data, sim.imp_scale.data)
                     mark_inactivelayer_bits!(sim.point_types.data)
                     verbose && println("Surface Refinement")
                     # Maximum spacing between (refined) surface ticks (if min_spacing = 1e-4m, only surface intervals wider than 0.1mm get new ticks)
-                    if sim.spacing_surface_refinement !== nothing
-                        refine_surface!(sim, sim.spacing_surface_refinement; update_other_fields=true)
-                    else
-                        refine_surface!(sim; update_other_fields=true)
-                    end
+                    refine_surface!(sim, sim.spacing_surface_refinement; update_other_fields=true)
+
                     update_till_convergence!( sim, potential_type, convergence_limit,
                                               n_iterations_between_checks = n_iterations_between_checks,
                                               max_n_iterations = max_n_iterations,

From 8d73f06992586fe7074c921b7156e6440a72cc35 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 24 Feb 2026 15:04:45 +0100
Subject: [PATCH 08/12] Add tests for edge cases

---
 test/test_refine_existing_potential.jl | 99 ++++++++++++++++++++++++++
 1 file changed, 99 insertions(+)

diff --git a/test/test_refine_existing_potential.jl b/test/test_refine_existing_potential.jl
index c65bd0ff6..aa4276566 100644
--- a/test/test_refine_existing_potential.jl
+++ b/test/test_refine_existing_potential.jl
@@ -120,3 +120,102 @@ end
     @test length(sim_cart_mm.electric_potential.grid[2]) == 60
     @test length(sim_cart_mm.electric_potential.grid[3]) == 64
 end
+
+@testset "Surface refinement edge cases" begin
+    T = Float64
+    # Surface model but no spacing_surface_refinement
+    begin
+        sim = @test_nowarn Simulation{T}(SSD_examples[:TrueCoaxial])
+        
+        @test_logs (:warn, r"Surface model detected but `spacing_surface_refinement` is not defined") begin
+            timed_calculate_electric_potential!(sim, verbose = false, depletion_handling = true)
+        end
+    end
+    
+    # Wrong refinement limits
+    begin
+        sim_cyl = @test_nowarn Simulation{T}(SSD_examples[:IVCIlayer])
+        sim_cyl.spacing_surface_refinement = (1e-3, 1e-3, 1e-3)
+        # Test normal behaviour Cylindrical
+        timed_calculate_electric_potential!(sim_cyl, verbose = false, depletion_handling = true)
+        grid_ax1_cyl = length(sim_cyl.electric_potential.grid[1])
+        grid_ax2_cyl = length(sim_cyl.electric_potential.grid[2])
+        grid_ax3_cyl = length(sim_cyl.electric_potential.grid[3])
+        
+        # --- Case: < 3 refinements ---
+        bad_limits = [0.2, 0.1]
+        @test_logs (:warn, r"Surface model detected") (:warn, r"Falling back to default") match_mode=:any begin
+            timed_calculate_electric_potential!(sim_cyl, refinement_limits=bad_limits, verbose = false, depletion_handling = true)
+        end
+        
+        # --- Case: last refinement > 0.05 ---
+        bad_limits = [0.2, 0.1, 0.09]
+        @test_logs (:warn, r"Surface model detected:") (:warn, r"Falling back to default") match_mode=:any begin
+            timed_calculate_electric_potential!(sim_cyl, refinement_limits=bad_limits, verbose = false, depletion_handling = true)
+        end
+
+        # Test normal behaviour Cartesian
+        config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
+        config_dict["grid"]["coordinates"] = "cartesian"
+        config_dict["grid"]["axes"]["x"] = Dict(
+            "from" => "-40",
+            "to"   => "40",
+            "boundaries"   => "inf"
+        )
+        config_dict["grid"]["axes"]["y"] = Dict(
+            "from" => "-40",
+            "to"   => "40",
+            "boundaries"   => "inf"
+        )
+        config_dict["grid"]["axes"]["z"] = Dict(
+            "from" => "-10",
+            "to"   => "90",
+            "boundaries"   => "inf"
+        )
+        sim_cart = @test_nowarn Simulation{T}(config_dict)
+        sim_cart.spacing_surface_refinement = (1e-3, 1e-3, 1e-3)
+        timed_calculate_electric_potential!(sim_cart, verbose = false, depletion_handling = true)
+        grid_ax1_cart = length(sim_cart.electric_potential.grid[1])
+        grid_ax2_cart = length(sim_cart.electric_potential.grid[2])
+	grid_ax3_cart = length(sim_cart.electric_potential.grid[3])
+        
+        @test grid_ax1_cyl == 74
+        @test grid_ax2_cyl == 1
+        @test grid_ax3_cyl == 143
+        @test grid_ax1_cart == 178
+        @test grid_ax2_cart == 178
+        @test grid_ax3_cart == 149
+    end
+    
+    # Spacing out of bounds
+    sim = @test_nowarn Simulation{T}(SSD_examples[:IVCIlayer])
+    large_spacing = (1e-2, 1e-2, 1e-2)
+    sim.spacing_surface_refinement = large_spacing
+    
+    @test_logs (
+        (:warn, r"out of bounds"s)
+    ) match_mode=:any begin
+        timed_calculate_electric_potential!(sim, verbose=false, depletion_handling=true)
+    end
+          
+    # Cylindrical grid: φ-axis not refined
+    begin
+        config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
+        
+        config_dict["grid"]["axes"]["phi"] = Dict(
+            "from" => "0",
+            "to"   => "120"
+        )
+        
+        sim_phi = @test_nowarn Simulation{T}(config_dict)
+        sim_phi.spacing_surface_refinement = nothing
+        timed_calculate_electric_potential!(sim_phi, verbose = false, depletion_handling = true)
+        notref_phi_len = length(sim_phi.electric_potential.grid.axes[2].ticks)
+
+        sim_phi.spacing_surface_refinement = (1e-3, 1e-3, 1e-3)
+        timed_calculate_electric_potential!(sim_phi, verbose = false, depletion_handling = true)
+        ref_phi_len = length(sim_phi.electric_potential.grid.axes[2].ticks)
+        
+        @test ref_phi_len == notref_phi_len
+    end
+end

From a0724bc44bc8f363b6f48684cd91d64e1fc6596e Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 12 Mar 2026 17:08:49 +0100
Subject: [PATCH 09/12] Pass surface refinement to World

---
 src/Simulation/Simulation.jl | 45 +++++++++++++-----------------------
 src/World/World.jl           | 31 +++++++++++++++----------
 2 files changed, 35 insertions(+), 41 deletions(-)

diff --git a/src/Simulation/Simulation.jl b/src/Simulation/Simulation.jl
index 000a5c6bf..1f5f9cd72 100644
--- a/src/Simulation/Simulation.jl
+++ b/src/Simulation/Simulation.jl
@@ -38,7 +38,6 @@ mutable struct Simulation{T <: SSDFloat, CS <: AbstractCoordinateSystem} <: Abst
     electric_potential::Union{ElectricPotential{T}, Missing}
     weighting_potentials::Vector{Any}
     electric_field::Union{ElectricField{T}, Missing}
-    spacing_surface_refinement::Union{NTuple{3,T}, Nothing}
 end
 
 function Simulation{T,CS}() where {T <: SSDFloat, CS <: AbstractCoordinateSystem}
@@ -55,8 +54,7 @@ function Simulation{T,CS}() where {T <: SSDFloat, CS <: AbstractCoordinateSystem
         missing,
         missing,
         [missing],
-        missing,
-        nothing
+        missing
     )
 end
 
@@ -115,7 +113,6 @@ function Simulation(nt::NamedTuple)
     else
         [missing for contact in sim.detector.contacts]
     end
-    sim.spacing_surface_refinement = get(nt, :spacing_surface_refinement, nothing)
     return sim
 end
 Base.convert(T::Type{Simulation}, x::NamedTuple) = T(x)
@@ -188,12 +185,6 @@ function Simulation{T}(dict::AbstractDict)::Simulation{T} where {T <: SSDFloat}
         end
     end
     sim.weighting_potentials = Missing[ missing for i in 1:length(sim.detector.contacts)]
-    sim.spacing_surface_refinement = if haskey(dict, "grid") && isa(dict["grid"], AbstractDict) &&
-                                       haskey(dict["grid"], "spacing_surface_refinement")
-        ntuple(i -> T(dict["grid"]["spacing_surface_refinement"][i]), 3)
-    else
-        nothing
-    end
     return sim
 end
 
@@ -778,8 +769,8 @@ function refine!(sim::Simulation{T}, ::Type{ElectricPotential},
     sim.electric_potential = refine_scalar_potential(sim.electric_potential, T.(max_diffs), T.(minimum_distances))
 
     if update_other_fields
-        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data,
-                                        not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts)
+        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data;
+                                        not_only_paint_contacts, paint_contacts)
 
         
         sim.imp_scale = ImpurityScale(ImpurityScaleArray(pcs), sim.electric_potential.grid)
@@ -824,7 +815,7 @@ function refine!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int
 end
 
 """
-refine_surface!(sim::Simulation{T};
+        refine_surface!(sim::Simulation{T};
                     min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4)),
                     not_only_paint_contacts::Bool = true,
                     paint_contacts::Bool = true,
@@ -846,13 +837,13 @@ Special behaviour:
 
 # Arguments
 - `sim::Simulation{T}`: The simulation object containing the electric potential to refine.
-- `min_spacing::NTuple{3,T}`: Minimum allowed spacing in the surface after refinement. Default is `(1e-4, 1e-4, 1e-4)m`.
+- `min_spacing::NTuple{3,T}`: Minimum allowed spacing in the surface after refinement.
 """
 
-function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4));
+function refine_surface!(sim::Simulation{T,CS}, min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4));
                          not_only_paint_contacts::Bool = true,
                          paint_contacts::Bool = true,
-                         update_other_fields::Bool = true) where {T <: SSDFloat}
+                         update_other_fields::Bool = true) where {T <: SSDFloat, CS <: AbstractCoordinateSystem}
     
     old_grid = sim.electric_potential.grid
 
@@ -888,7 +879,7 @@ function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4)
     end, 3)
 
     # Refine axes: skip phi-axis if Cylindrical
-    is_cyl = old_grid isa Grid{<:Any, 3, Cylindrical}
+    is_cyl = CS === Cylindrical
     new_axes = ntuple(i -> begin
                           if is_cyl && i == 2
                               # Keep φ-axis type but close it for interpolation
@@ -897,13 +888,9 @@ function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4)
                               _refine_axis_surface( old_grid.axes[i], surface_intervals[i], min_spacing[i],)
                           end
                       end, 3)
-    
-    if is_cyl
-        new_grid = Grid{T,3,Cylindrical}((new_axes[1], new_axes[2], new_axes[3]))
-    else
-        new_grid = Grid{T,3,Cartesian}((new_axes[1], new_axes[2], new_axes[3]))
-    end
 
+    new_grid = Grid{T, 3, CS}((new_axes[1], new_axes[2], new_axes[3]))
+    
     # Interpolate potential onto new grid
     closed_pot = _get_closed_potential(sim.electric_potential)
     new_data = Array{T,3}(undef, size(new_grid))
@@ -930,8 +917,8 @@ function refine_surface!(sim::Simulation{T}, min_spacing::NTuple{3,T} = (T(1e-4)
     # Update dependent fields
     if update_other_fields
 
-        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data,
-                                        not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts)
+        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data;
+                                        not_only_paint_contacts, paint_contacts)
         sim.imp_scale = ImpurityScale(ImpurityScaleArray(pcs), sim.electric_potential.grid)
         sim.q_eff_imp = EffectiveChargeDensity(EffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
         sim.q_eff_fix = EffectiveChargeDensity(FixedEffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
@@ -1033,8 +1020,8 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
         end
         has_surface_model = isdefined(sim.detector.semiconductor.impurity_density_model, :surface_imp_model)
 
-        if has_surface_model    
-            if sim.spacing_surface_refinement === nothing
+        if has_surface_model
+            if isnothing(sim.world.spacing_surface_refinement)
                 @warn """Surface model detected but `spacing_surface_refinement` is not defined.
                          Surface refinement will not be performed."""
             else
@@ -1171,13 +1158,13 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                                                 paint_contacts = paint_contacts,
                                                 sor_consts = is_last_ref ? T(1) : sor_consts )
                 
-                if has_surface_model && iref==3 && sim.spacing_surface_refinement !== nothing
+                if has_surface_model && iref==3 && !isnothing(sim.world.spacing_surface_refinement)
                     mark_bulk_bits!(sim.point_types.data)
                     mark_undep_bits!(sim.point_types.data, sim.imp_scale.data)
                     mark_inactivelayer_bits!(sim.point_types.data)
                     verbose && println("Surface Refinement")
                     # Maximum spacing between (refined) surface ticks (if min_spacing = 1e-4m, only surface intervals wider than 0.1mm get new ticks)
-                    refine_surface!(sim, sim.spacing_surface_refinement; update_other_fields=true)
+                    refine_surface!(sim, sim.world.spacing_surface_refinement; update_other_fields=true)
 
                     update_till_convergence!( sim, potential_type, convergence_limit,
                                               n_iterations_between_checks = n_iterations_between_checks,
diff --git a/src/World/World.jl b/src/World/World.jl
index 093c15c9d..1b8299316 100644
--- a/src/World/World.jl
+++ b/src/World/World.jl
@@ -55,18 +55,25 @@ Definition of the finite volume on which a [`Simulation`](@ref) is performed.
 """
 struct World{T <: SSDFloat, N, S} <: AbstractWorld{T, N} 
     intervals::NTuple{N, SSDInterval{T}}
+    spacing_surface_refinement::Union{NTuple{N,T}, Nothing}
 end
 
 function World{T, N, S}(args...) where {T <: SSDFloat, N, S} 
     return World{T, N, S}(args)
 end
 
-
 function World(T, dict::AbstractDict, input_units::NamedTuple)::World
+
+    spacing_surface_refinement = if haskey(dict, "spacing_surface_refinement")
+        ntuple(i -> T(dict["spacing_surface_refinement"][i]), 3)
+    else
+        nothing
+    end
+
     if dict["coordinates"] == "cylindrical"
-        CylindricalWorld(T, dict["axes"], input_units)
+        return CylindricalWorld(T, dict["axes"], input_units, spacing_surface_refinement)
     elseif dict["coordinates"] == "cartesian"
-        CartesianWorld(T, dict["axes"], input_units)
+        return CartesianWorld(T, dict["axes"], input_units, spacing_surface_refinement)
     else
         error("Gridtype must be \"cylindrical\" or \"cartesian\"")
     end
@@ -158,37 +165,37 @@ function get_cartesian_SSDInterval(T, dict::AbstractDict, input_units::NamedTupl
 end
 
 
-function CylindricalWorld(T, dict::AbstractDict, input_units::NamedTuple)::World
+function CylindricalWorld(T, dict::AbstractDict, input_units::NamedTuple, spacing_surface_refinement)::World
     r_int = get_r_SSDInterval(T, dict["r"], input_units)
     φ_int = get_φ_SSDInterval(T, dict, input_units)
     z_int = get_cartesian_SSDInterval(T, dict["z"], input_units)
-    return World{T, 3, Cylindrical}( r_int, φ_int, z_int )
+    return World{T, 3, Cylindrical}( (r_int, φ_int, z_int), spacing_surface_refinement )
 end
 
 
-function CartesianWorld(T, dict::AbstractDict, input_units::NamedTuple)::World
+function CartesianWorld(T, dict::AbstractDict, input_units::NamedTuple, spacing_surface_refinement)::World
     x_int = get_cartesian_SSDInterval(T, dict["x"], input_units)
     y_int = get_cartesian_SSDInterval(T, dict["y"], input_units)
     z_int = get_cartesian_SSDInterval(T, dict["z"], input_units)
-    return World{T, 3, Cartesian}( x_int, y_int, z_int )
+    return World{T, 3, Cartesian}( (x_int, y_int, z_int), spacing_surface_refinement )
 end
 
-function CartesianWorld(xl::T, xr::T, yl::T, yr::T, zl::T, zr::T)::World where {T <: SSDFloat}
+function CartesianWorld(xl::T, xr::T, yl::T, yr::T, zl::T, zr::T; spacing_surface_refinement=nothing)::World where {T <: SSDFloat}
     Δx::T = (xr - xl) / 10
     Δy::T = (yr - yl) / 10
     Δz::T = (zr - zl) / 10
     x_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(xl - Δx, xr + Δx)
     y_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(yl - Δy, yr + Δy)
     z_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(zl - Δz, zr + Δz)
-    return World{T, 3, Cartesian}( x_int, y_int, z_int )
+    return World{T, 3, Cartesian}( (x_int, y_int, z_int), spacing_surface_refinement )
 end
 
-function CylindricalWorld(r_max::T, zl::T, zr::T)::World where {T <: SSDFloat}
+function CylindricalWorld(r_max::T, zl::T, zr::T; spacing_surface_refinement=nothing)::World where {T <: SSDFloat}
     r_int = SSDInterval{T, :closed, :closed, :r0, :infinite}(T(0), abs(r_max * T(1.1)))
     φ_int = SSDInterval{T, :closed, :open, :periodic, :periodic}(T(0), T(2π))
     Δz::T = zr - zl
     z_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(zl - (Δz / 10), zr + (Δz / 10))
-    return World{T, 3, Cylindrical}( r_int, φ_int, z_int )
+    return World{T, 3, Cylindrical}( (r_int, φ_int, z_int), spacing_surface_refinement )
 end
 
 function World(::Type{Cylindrical}, limits::NTuple{6, T})::World where {T <: SSDFloat}
@@ -205,4 +212,4 @@ end
 function max_tick_distance_default(w::World{T, 3, Cartesian}; n = 50) where {T}
     Δw = width.(w.intervals)
     return max(Δw...) * internal_length_unit / n
-end
\ No newline at end of file
+end

From 4fc9950271475a41b1b55cd5f1772cfb9b5ee7c2 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 12 Mar 2026 17:09:20 +0100
Subject: [PATCH 10/12] suggested change

---
 src/PotentialCalculation/Refinement.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/PotentialCalculation/Refinement.jl b/src/PotentialCalculation/Refinement.jl
index 055ed0a7f..9ef1a94de 100644
--- a/src/PotentialCalculation/Refinement.jl
+++ b/src/PotentialCalculation/Refinement.jl
@@ -191,9 +191,7 @@ function _refine_axis_surface( ax::DiscreteAxis{T, <:Any, <:Any, ClosedInterval{
     first_surface = findfirst(surface_intervals)
     last_surface  = findlast(surface_intervals)
 
-    if first_surface === nothing || last_surface === nothing
-        return ax # Nothing to refine
-    end
+    (isnothing(first_surface) || isnothing(last_surface)) && return ax
 
     # Create flags for all intervals to refine
     refine_flags = copy(surface_intervals)

From 657f72fc20aab5a9ca5d33289edc56e8f8fd8512 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 12 Mar 2026 17:09:34 +0100
Subject: [PATCH 11/12] Change tests to new format

---
 test/test_refine_existing_potential.jl | 31 ++++++++++++++------------
 1 file changed, 17 insertions(+), 14 deletions(-)

diff --git a/test/test_refine_existing_potential.jl b/test/test_refine_existing_potential.jl
index aa4276566..2dd34ce36 100644
--- a/test/test_refine_existing_potential.jl
+++ b/test/test_refine_existing_potential.jl
@@ -134,8 +134,10 @@ end
     
     # Wrong refinement limits
     begin
-        sim_cyl = @test_nowarn Simulation{T}(SSD_examples[:IVCIlayer])
-        sim_cyl.spacing_surface_refinement = (1e-3, 1e-3, 1e-3)
+        config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
+        config_dict["grid"]["spacing_surface_refinement"] = [1e-3, 1e-3, 1e-3]
+        sim_cyl = @test_nowarn Simulation{T}(config_dict)
+
         # Test normal behaviour Cylindrical
         timed_calculate_electric_potential!(sim_cyl, verbose = false, depletion_handling = true)
         grid_ax1_cyl = length(sim_cyl.electric_potential.grid[1])
@@ -155,7 +157,6 @@ end
         end
 
         # Test normal behaviour Cartesian
-        config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
         config_dict["grid"]["coordinates"] = "cartesian"
         config_dict["grid"]["axes"]["x"] = Dict(
             "from" => "-40",
@@ -173,7 +174,6 @@ end
             "boundaries"   => "inf"
         )
         sim_cart = @test_nowarn Simulation{T}(config_dict)
-        sim_cart.spacing_surface_refinement = (1e-3, 1e-3, 1e-3)
         timed_calculate_electric_potential!(sim_cart, verbose = false, depletion_handling = true)
         grid_ax1_cart = length(sim_cart.electric_potential.grid[1])
         grid_ax2_cart = length(sim_cart.electric_potential.grid[2])
@@ -188,9 +188,9 @@ end
     end
     
     # Spacing out of bounds
-    sim = @test_nowarn Simulation{T}(SSD_examples[:IVCIlayer])
-    large_spacing = (1e-2, 1e-2, 1e-2)
-    sim.spacing_surface_refinement = large_spacing
+    config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
+    config_dict["grid"]["spacing_surface_refinement"] = [1e-2, 1e-2, 1e-2]
+    sim = @test_nowarn Simulation{T}(config_dict)
     
     @test_logs (
         (:warn, r"out of bounds"s)
@@ -208,13 +208,16 @@ end
         )
         
         sim_phi = @test_nowarn Simulation{T}(config_dict)
-        sim_phi.spacing_surface_refinement = nothing
-        timed_calculate_electric_potential!(sim_phi, verbose = false, depletion_handling = true)
-        notref_phi_len = length(sim_phi.electric_potential.grid.axes[2].ticks)
-
-        sim_phi.spacing_surface_refinement = (1e-3, 1e-3, 1e-3)
-        timed_calculate_electric_potential!(sim_phi, verbose = false, depletion_handling = true)
-        ref_phi_len = length(sim_phi.electric_potential.grid.axes[2].ticks)
+        
+        delete!(config_dict["grid"], "spacing_surface_refinement")
+        sim_phi_noref = @test_nowarn Simulation{T}(config_dict)
+        timed_calculate_electric_potential!(sim_phi_noref, verbose = false, depletion_handling = true)
+        notref_phi_len = length(sim_phi_noref.electric_potential.grid.axes[2].ticks)
+
+        config_dict["grid"]["spacing_surface_refinement"] = [1e-3,1e-3,1e-3]
+        sim_phi_ref = @test_nowarn Simulation{T}(config_dict)
+        timed_calculate_electric_potential!(sim_phi_ref, verbose = false, depletion_handling = true)
+        ref_phi_len = length(sim_phi_ref.electric_potential.grid.axes[2].ticks)
         
         @test ref_phi_len == notref_phi_len
     end

From 23c3833d4a4547e9857cb821a08d34566edfafd2 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Fri, 13 Mar 2026 10:30:30 +0100
Subject: [PATCH 12/12] adding is_cyl and _parse_value

---
 src/Simulation/Simulation.jl | 2 +-
 src/World/World.jl           | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/Simulation/Simulation.jl b/src/Simulation/Simulation.jl
index 1f5f9cd72..161bc3906 100644
--- a/src/Simulation/Simulation.jl
+++ b/src/Simulation/Simulation.jl
@@ -894,7 +894,7 @@ function refine_surface!(sim::Simulation{T,CS}, min_spacing::NTuple{3,T} = (T(1e
     # Interpolate potential onto new grid
     closed_pot = _get_closed_potential(sim.electric_potential)
     new_data = Array{T,3}(undef, size(new_grid))
-    only2d = size(closed_pot.data, 2) == 1
+    only2d = is_cyl && size(closed_pot.data, 2) == 1
     int = interpolate_closed_potential(closed_pot, Val(only2d))
     
     if only2d
diff --git a/src/World/World.jl b/src/World/World.jl
index 1b8299316..b382b0d24 100644
--- a/src/World/World.jl
+++ b/src/World/World.jl
@@ -63,13 +63,13 @@ function World{T, N, S}(args...) where {T <: SSDFloat, N, S}
 end
 
 function World(T, dict::AbstractDict, input_units::NamedTuple)::World
-
+    
     spacing_surface_refinement = if haskey(dict, "spacing_surface_refinement")
-        ntuple(i -> T(dict["spacing_surface_refinement"][i]), 3)
+        ntuple(i -> _parse_value(T, dict["spacing_surface_refinement"][i], internal_length_unit), 3)
     else
         nothing
     end
-
+    
     if dict["coordinates"] == "cylindrical"
         return CylindricalWorld(T, dict["axes"], input_units, spacing_surface_refinement)
     elseif dict["coordinates"] == "cartesian"