From 2024a8bede5bda1ae17d9c40ece2453316b95e0b Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Mon, 31 Aug 2026 12:00:30 -0400 Subject: [PATCH] VACUUM - IMPROVEMENT - adding 3D conformal walls --- src/Vacuum/DataTypes.jl | 107 +++++++++++++++++++++++++++++++------ src/Vacuum/Vacuum.jl | 2 +- test/runtests_vacuum.jl | 113 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 203 insertions(+), 19 deletions(-) diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index da0b8167f..de1acfe02 100644 --- a/src/Vacuum/DataTypes.jl +++ b/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, diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 3fd13eb22..9c4921618 100644 --- a/src/Vacuum/Vacuum.jl +++ b/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 diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 19f80b94d..2db363573 100644 --- a/test/runtests_vacuum.jl +++ b/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")