Vacuum - FEATURE! - Add the 3D surface-current matrix and halve the solve with stellarator symmetry

src/Vacuum/DataTypes.jl

@@ -154,7 +154,7 @@ boundary-integral solve produces along the way.
 
   - `wv::Matrix{ComplexF64}`: Vacuum energy matrix Wᵛ (`num_modes × num_modes`), block-diagonal in n for 2D
   - `I_v::Matrix{ComplexF64}`: Vacuum surface-current matrix Iᵛ (`num_modes × num_modes`), left zeroed
-    unless `compute_vacuum_response` is called with `compute_Iv=true` (2D only). Stored without the
+    unless `compute_vacuum_response` is called with `compute_Iv=true`. Stored without the
     `μ₀`/`4π²` normalization: the physical surface inductance is `μ₀(2π)²·I_v⁻¹`
     (see `PerturbedEquilibrium.calc_surface_inductance`).
   - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`)
@@ -198,6 +198,9 @@ Struct containing input settings for vacuum wall geometry.
         + `"mod_dee"`: Modified Dee-shaped wall
         + `"filepath"`: Custom wall shape from the file you specify
 
+    A non-axisymmetric boundary (`nzeta_in > 1`) supports `"nowall"` and `"conformal"`; the others
+    are poloidal contours that get revolved and so need `nzeta_in == 1`.
+
   - `a::Float64`: Distance of wall from plasma in units of major radius (conformal), or shape parameter (others)
 
   - `aw::Float64`: Half-thickness parameter for Dee-shaped walls
@@ -213,7 +216,9 @@ Struct containing input settings for vacuum wall geometry.
     # Core shape selection
 
   - `equal_arc_wall::Bool`: Flag to enforce equal arc length distribution of nodes on the wall
-    (recommended unless wall is very close to plasma)
+    (recommended unless wall is very close to plasma). Re-parameterizing the contour breaks the
+    plasma/wall grid index correspondence that the singular quadrature relies on; it is ignored
+    outright for `nzeta_in > 1` walls.
 """
 @kwdef struct WallShapeSettings
 
@@ -600,8 +605,8 @@ end
 """
     WallGeometry3D
 
-Struct holding wall geometry data for vacuum calculations. Arrays are of length
-`mtheta`, where `mtheta` is the number of poloidal grid points and θ ∈ [0, 1).
+Struct holding wall geometry data for vacuum calculations. Arrays have `mtheta * nzeta` rows,
+ordered with the poloidal index fastest.
 
 # Fields
 
@@ -627,21 +632,36 @@ struct WallGeometry3D
 end
 
 """
-    WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings)
+    WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall_settings::WallShapeSettings) -> WallGeometry3D
 
-Constructor to initialize the 3D wall geometry based on the provided vacuum inputs and wall shape settings.
-Currently only works for axisymmetric walls generated by toroidal extrusion of 2D poloidal contours.
+Constructor to initialize the 3D wall geometry based on the provided vacuum inputs, plasma surface and
+wall shape settings. This is the 3D counterpart of [`WallGeometry`](@ref) and selects the shape the same
+way: an axisymmetric boundary builds the 2D poloidal contour and revolves it, so every shape
+[`WallGeometry`](@ref) offers is available; a non-axisymmetric boundary offsets the plasma surface.
+
+Expects a full-torus boundary — call [`expand_field_periods`](@ref) first when `nfp > 1`.
 
 # Arguments
 
   - `inputs::VacuumInput`: Struct containing vacuum calculation parameters
+  - `plasma_surf::PlasmaGeometry3D`: Plasma surface the wall is built around (used by the conformal wall)
   - `wall_settings::WallShapeSettings`: Struct specifying wall shape and parameters
 
 # Returns
 
   - `WallGeometry3D`: Struct containing wall surface coordinates and derivatives
+
+# Notes
+
+  - Axisymmetric boundaries (`nzeta_in == 1`) support nowall, conformal, elliptical, dee, mod_dee and
+    from_file; non-axisymmetric boundaries support nowall and conformal
+  - The conformal wall for a non-axisymmetric boundary displaces each plasma point along its own normal,
+    which keeps wall grid index `(i, j)` the closest wall point to plasma index `(i, j)` — the
+    correspondence the near-field patch of `compute_3D_kernel_matrices!` assumes
+  - Rejects a conformal offset that folds the surface, and warns when the plasma–wall gap is smaller
+    than one cell of the coarser grid (the double-layer near-field quadrature cannot resolve it)
 """
-function WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings)
+function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall_settings::WallShapeSettings)
 
     # Basic wall flags
     nowall = wall_settings.shape == "nowall"
@@ -676,16 +696,35 @@ function WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings)
         )
     end
 
-    inputs.nzeta_in > 1 && error("3D wall geometry not yet implemented for non-axisymmetric inputs")
-    inputs.nfp > 1 && error("3D wall geometry not yet implemented for nfp > 1")
-
-    # Plasma surface coordinates (2D)
-    surf_2D = PlasmaGeometry(inputs)
-    wall_2D = WallGeometry(inputs, surf_2D, wall_settings)
+    # expand_field_periods leaves an axisymmetric boundary alone, so nfp > 1 can only arrive here with
+    # per-period counts that the rest of the 3D solve reads as full-torus ones.
+    inputs.nfp > 1 && error("3D wall geometry requires a full-torus boundary. Call expand_field_periods(inputs) first, or supply the whole torus with nfp = 1.")
 
-    # Build 3D surface point-by-point from 2D contour
-    for i in 1:mtheta, (j, ϕ) in enumerate(ϕ_grid)
-        r[i+mtheta*(j-1), :] .= [wall_2D.x[i] * cos(ϕ), wall_2D.x[i] * sin(ϕ), wall_2D.z[i]]
+    if inputs.nzeta_in == 1
+        # Axisymmetric boundary: build the 2D poloidal contour and revolve it toroidally
+        surf_2D = PlasmaGeometry(inputs)
+        wall_2D = WallGeometry(inputs, surf_2D, wall_settings)
+        for i in 1:mtheta, (j, ϕ) in enumerate(ϕ_grid)
+            r[i+mtheta*(j-1), :] .= [wall_2D.x[i] * cos(ϕ), wall_2D.x[i] * sin(ϕ), wall_2D.z[i]]
+        end
+    elseif wall_settings.shape == "conformal"
+        # Displace every plasma point outward along its own normal
+        wall_settings.equal_arc_wall && @warn "equal_arc_wall is ignored for non-axisymmetric (nzeta_in > 1) walls: it re-parameterizes a 2D contour and would break the plasma/wall index alignment the near-field patch relies on."
+
+        # Same logic as 2D
+        R_plasma = [hypot(plasma_surf.r[idx, 1], plasma_surf.r[idx, 2]) for idx in axes(plasma_surf.r, 1)]
+        offset_gap = wall_settings.a * 0.5 * (maximum(R_plasma) - minimum(R_plasma))
+        @info "Calculating conformal wall shape $((@sprintf "%.2e" offset_gap)) m from plasma surface."
+
+        # Plasma normal points into the plasma, so the displacement is along -normal
+        @inbounds for idx in axes(plasma_surf.r, 1)
+            scale = -offset_gap / sqrt(plasma_surf.normal[idx, 1]^2 + plasma_surf.normal[idx, 2]^2 + plasma_surf.normal[idx, 3]^2)
+            for k in 1:3
+                r[idx, k] = plasma_surf.r[idx, k] + scale * plasma_surf.normal[idx, k]
+            end
+        end
+    else
+        error("Wall shape $(wall_settings.shape) is not available for a non-axisymmetric boundary (nzeta_in > 1).")
     end
 
     # Compute tangent vectors and normal vectors via periodic bicubic splines
@@ -705,6 +744,42 @@ function WallGeometry3D(inputs::VacuumInput, wall_settings::WallShapeSettings)
     normal_orient = normal[idx, 1] > 0 ? 1 : -1
     @views normal .*= normal_orient
 
+    # Same-index plasma–wall separation: the near-field patch is centred on this pair.
+    gap = minimum(hypot(r[idx, 1] - plasma_surf.r[idx, 1], r[idx, 2] - plasma_surf.r[idx, 2], r[idx, 3] - plasma_surf.r[idx, 3]) for idx in axes(r, 1))
+
+    # Fold check needs a pointwise normal offset (same-index pair) - equal_arc_wall re-parameterizes and breaks that pairing
+    if wall_settings.shape == "conformal" && (inputs.nzeta_in > 1 || !wall_settings.equal_arc_wall)
+        # Offset stays regular while wall and plasma normals stay aligned and the area element has not collapsed
+        # Both fail when the gap exceeds the local concave radius of curvature.
+        min_align = Inf # min n̂_wall · n̂_plasma_outward; +1 healthy, ≤0 folded
+        min_area_ratio = Inf # min ||n_wall||/||n_plasma||; →0 at a caustic
+        @inbounds for idx in axes(normal, 1)
+            norm_wall = sqrt(normal[idx, 1]^2 + normal[idx, 2]^2 + normal[idx, 3]^2)
+            norm_plasma = sqrt(plasma_surf.normal[idx, 1]^2 + plasma_surf.normal[idx, 2]^2 + plasma_surf.normal[idx, 3]^2)
+            # Plasma n in, wall n out of vacuum → healthy offset dots to −1; store −cosine so +1 is aligned.
+            dot_wp = normal[idx, 1] * plasma_surf.normal[idx, 1] + normal[idx, 2] * plasma_surf.normal[idx, 2] + normal[idx, 3] * plasma_surf.normal[idx, 3]
+            min_align = min(min_align, -dot_wp / (norm_wall * norm_plasma))
+            min_area_ratio = min(min_area_ratio, norm_wall / norm_plasma)
+        end
+
+        min_align <= 0.0 && error(
+            "Conformal wall offset of $((@sprintf "%.2e" gap)) m self-intersects: it exceeds the local concave radius " *
+            "of curvature of the plasma surface somewhere, folding the offset surface. Reduce the wall distance a."
+        )
+        (min_align < 0.5 || min_area_ratio < 0.5) &&
+            @warn "Conformal wall is close to folding: min(n̂_wall·n̂_plasma)=$((@sprintf "%.3f" min_align)), min area-element ratio=$((@sprintf "%.3f" min_area_ratio)). Reduce the wall distance a."
+    end
+
+    # The 3D singular quadrature resolves the double-layer kernel well only while the plasma-wall gap spans a grid cell or more
+    h_θ = max(sqrt(sum(abs2, dr_dθ) / num_points), sqrt(sum(abs2, plasma_surf.dr_dθ) / num_points)) * dθ
+    h_ζ = max(sqrt(sum(abs2, dr_dζ) / num_points), sqrt(sum(abs2, plasma_surf.dr_dζ) / num_points)) * dζ
+    cells = gap / max(h_θ, h_ζ)
+    if cells < 1.0
+        @warn "Plasma–wall gap $((@sprintf "%.2e" gap)) m spans only $((@sprintf "%.2f" cells)) cells of the coarser grid " *
+                "(h_θ=$((@sprintf "%.2e" h_θ)) m, h_ζ=$((@sprintf "%.2e" h_ζ)) m). The double-layer near-field quadrature is inaccurate " *
+                "below one cell. Refine the coarser of mtheta/nzeta, or move the wall out."
+    end
+
     return WallGeometry3D(
         nowall,
         is_closed_toroidal,