Vacuum - FEATURE - Add 3D conformal conducting wall capabilities

src/Vacuum/DataTypes.jl

@@ -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,

src/Vacuum/Vacuum.jl

@@ -192,7 +192,7 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the
     # Full-torus geometry for source surface; observers are restricted to one field period
     full = expand_field_periods(inputs)
     plasma_surf = PlasmaGeometry3D(full)
-    wall = WallGeometry3D(full, wall_settings)
+    wall = WallGeometry3D(full, plasma_surf, wall_settings)
 
     num_points_per_fp = mtheta * nzeta # points per field period
     num_points = num_points_per_fp * nfp # full-torus point count

test/runtests_vacuum.jl

@@ -594,7 +594,7 @@
 
         @testset "WallGeometry3D nowall" begin
             inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17)
-            wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, WallShapeSettings(shape="nowall"))
+            wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs), WallShapeSettings(shape="nowall"))
             @test wall.nowall == true
             @test wall.mtheta == 32
             @test wall.nzeta == 32
@@ -603,7 +603,11 @@
 
         @testset "WallGeometry3D conformal" begin
             inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17)
-            wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, WallShapeSettings(shape="conformal", a=0.2))
+            wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(
+                inputs,
+                GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs),
+                WallShapeSettings(shape="conformal", a=0.2)
+            )
             @test wall.nowall == false
             @test wall.mtheta == 32
             @test wall.nzeta == 32
@@ -614,6 +618,111 @@
             @test all(isfinite, wall.normal)
         end
 
+        # Corrugated torus on a genuinely periodic (endpoint-excluded) grid. The shared
+        # _make_3d_nonaxis_inputs helper samples range(0, 2π, length=n), which repeats the seam
+        # point and leaves the surface non-smooth there — harmless for a nowall response, but it
+        # makes the offset surface fold, so the wall tests build their own boundary.
+        _make_3d_periodic_inputs(; mtheta=24, nzeta=24, mtheta_in=16, nzeta_in=16) = begin
+            θ_in = range(; start=0, length=mtheta_in, step=2π/mtheta_in)
+            ζ_in = range(; start=0, length=nzeta_in, step=2π/nzeta_in)
+            X = zeros(mtheta_in, nzeta_in)
+            Y = similar(X)
+            Z = similar(X)
+            for (i, θ) in enumerate(θ_in), (j, ζ) in enumerate(ζ_in)
+                R = 1.7 + 0.3 * cos(θ) + 0.05 * cos(2ζ) * cos(θ)
+                X[i, j] = R * cos(ζ)
+                Y[i, j] = R * sin(ζ)
+                Z[i, j] = 0.3 * sin(θ) + 0.05 * sin(2ζ) * sin(θ)
+            end
+            VacuumInput(x=vec(X), y=vec(Y), z=vec(Z), mtheta_in=mtheta_in, nzeta_in=nzeta_in,
+                m_modes=[1, 2], n_modes=[0, 1], mtheta=mtheta, nzeta=nzeta)
+        end
+
+        @testset "WallGeometry3D conformal, non-axisymmetric input" begin
+            inputs = _make_3d_periodic_inputs(mtheta=24, nzeta=24)
+            settings = WallShapeSettings(shape="conformal", a=0.2, equal_arc_wall=false)
+            plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs)
+            wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, settings)
+
+            num_points = 24 * 24
+            @test wall.nowall == false
+            @test size(wall.r) == (num_points, 3)
+            @test all(isfinite, wall.r) && all(isfinite, wall.normal)
+
+            # Uniform offset along the plasma normal: every point moves the same distance, outward.
+            offsets = [norm(wall.r[i, :] - plasma.r[i, :]) for i in 1:num_points]
+            @test maximum(offsets) - minimum(offsets) < 1e-12
+            @test all(hypot.(wall.r[:, 1], wall.r[:, 2]) .> 0)
+            R_wall = [hypot(wall.r[i, 1], wall.r[i, 2]) for i in 1:num_points]
+            R_plasma = [hypot(plasma.r[i, 1], plasma.r[i, 2]) for i in 1:num_points]
+            @test maximum(R_wall) > maximum(R_plasma)
+
+            # The offset point is the closest wall point to its own plasma point, which is the
+            # index alignment the near-field patch assumes.
+            for i in (1, 300, num_points)
+                @test offsets[i] ≈ minimum(norm(wall.r[j, :] - plasma.r[i, :]) for j in 1:num_points)
+            end
+
+            # Wall normals face out of the vacuum region, opposite the inward plasma normals.
+            aligns = [dot(wall.normal[i, :], plasma.normal[i, :]) for i in 1:num_points]
+            @test all(aligns .< 0)
+        end
+
+        @testset "WallGeometry3D non-axisymmetric error paths" begin
+            inputs = _make_3d_periodic_inputs(mtheta=24, nzeta=24)
+            plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs)
+            for shape in ("elliptical", "dee", "mod_dee", "some_wall_file.dat")
+                @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, WallShapeSettings(shape=shape))
+            end
+            # equal_arc_wall has no meaning without a 2D contour and is ignored with a warning
+            @test_logs (:warn, r"equal_arc_wall is ignored") match_mode=:any GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(
+                inputs,
+                plasma,
+                WallShapeSettings(shape="conformal", a=0.2, equal_arc_wall=true)
+            )
+            # a full-torus boundary is required; expand_field_periods must run first
+            per_period = VacuumInput(x=inputs.x, y=inputs.y, z=inputs.z, mtheta_in=16, nzeta_in=16,
+                m_modes=inputs.m_modes, n_modes=inputs.n_modes, mtheta=24, nzeta=24, nfp=3)
+            @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(per_period, plasma, WallShapeSettings(shape="conformal", a=0.2))
+            # an offset that folds the surface is rejected rather than silently returned
+            folded = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12)
+            @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(
+                folded,
+                GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(folded),
+                WallShapeSettings(shape="conformal", a=0.2, equal_arc_wall=false)
+            )
+        end
+
+        @testset "WallGeometry3D shape selection" begin
+            # Every shape reaches the right builder through the same chain WallGeometry uses in 2D.
+            # elongated so the elliptical branch has a well-defined focal distance
+            θ_eq = range(0, 2π, length=33)[1:32]
+            axi = VacuumInput(mtheta_in=32, nzeta_in=1, x=collect(1.7 .+ 0.3 .* cos.(θ_eq)),
+                z=collect(0.45 .* sin.(θ_eq)), ν=zeros(32), m_modes=[1, 2], n_modes=[0, 1],
+                mtheta=32, nzeta=32)
+            for (shape, settings) in (("elliptical", WallShapeSettings(shape="elliptical", a=0.5)),
+                ("dee", WallShapeSettings(shape="dee", a=0.3)),
+                ("mod_dee", WallShapeSettings(shape="mod_dee", a=0.5, cw=1.7)))
+                wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(axi, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(axi), settings)
+                @test size(wall.r) == (32 * 32, 3)
+                @test all(isfinite, wall.normal)
+                # revolved shapes are axisymmetric: R and Z repeat from one toroidal plane to the next
+                @test hypot(wall.r[1, 1], wall.r[1, 2]) ≈ hypot(wall.r[1+32, 1], wall.r[1+32, 2])
+                @test wall.r[1, 3] ≈ wall.r[1+32, 3]
+            end
+            # Gap warning uses the measured same-index separation, not a conformal-only offset that
+            # is 0 for these shapes. a = 1 m sits well outside one coarser cell, so no sub-cell warning.
+            @test_logs min_level=Base.CoreLogging.Warn GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(
+                axi, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(axi), WallShapeSettings(shape="elliptical", a=1.0)
+            )
+            # an unreadable wall file is still reported by the 2D reader
+            @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(
+                axi,
+                GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(axi),
+                WallShapeSettings(shape="no_such_wall.dat")
+            )
+        end
+
         @testset "compute_vacuum_response 3D nowall" begin
             inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17)
             wall_settings = WallShapeSettings(shape="nowall")