Utilities - FEATURE - Write the response, tearing growth rates and NTV torque to IMAS

src/GeneralizedPerturbedEquilibrium.jl

@@ -270,8 +270,21 @@ end
 """
     write_imas(dd, result)
 
-Write one `dd.mhd_linear` `toroidal_mode` entry per toroidal mode number, with the least-stable
-(minimum real part) free-boundary `energy_perturbed` of that n. `result` is the named tuple `main` returns.
+Write the key results of a run into the IMAS data dictionary `dd`. `result` is the named tuple `main` returns;
+a stage that did not run writes nothing.
+
+  - `dd.mhd_linear` `toroidal_mode[j]`, one per toroidal mode number n: the least-stable (minimum real part)
+    free-boundary `energy_perturbed` of that n and, after a perturbed-equilibrium solve, the response's normal
+    displacement ξ_n [m] (`plasma.displacement_perpendicular`) and normal field b_n [T]
+    (`plasma.b_field_perturbed.coordinate1`) by poloidal harmonic on the code's ψ grid
+    (`plasma.grid.dim1` = √ψ_N, `dim2` = m, grid type 24).
+  - Further `toroidal_mode` entries after those, `perturbation_type.name = "tearing"`, one per tearing root
+    (one per surface uncoupled, one coupled): `growthrate` [1/s] and `frequency` [Hz].
+  - `dd.core_sources` `source`, identifier `neoclassical`, one per kinetic-torque method: the NTV torque density
+    `momentum_tor` [N m⁻²] and the torque inside each surface `torque_tor_inside` [N m] on `grid.psi_norm`, with
+    `grid.rho_tor_norm` interpolated from `dd.equilibrium` when the dd holds the equilibrium.
+
+Δ' and the resonant fields at the rational surfaces have no IMAS slot and are not written.
 """
 function write_imas(dd, result)
     ffs = result.ffs
@@ -298,9 +311,83 @@ function write_imas(dd, result)
         mode.energy_perturbed = minimum(real.(free_energies.et[n_indices])) # least-stable energy for this n-toroidal mode
     end
 
+    pe = get(result, :pe, nothing)
+    pe === nothing || _write_imas_response!(ts, ffs, pe)
+    tearing = get(result, :tearing, nothing)
+    tearing === nothing || _write_imas_tearing!(ts, ffs, tearing)
+    kinetic = get(result, :kinetic, nothing)
+    kinetic === nothing || _write_imas_torque!(dd, result.equil, kinetic)
     return dd
 end
 
+# Response ξ_n and b_n of each n by poloidal harmonic, on √ψ_N × m (straight-field-line Fourier grid, type 24).
+function _write_imas_response!(ts, ffs, pe)
+    (pe.xi_n_modes === nothing || pe.b_n_modes === nothing || isempty(pe.psi_grid)) && return nothing
+    ms = collect(Float64, ffs.mlow:ffs.mhigh)
+    for j in 0:(ffs.npert-1)
+        cols = (j*ffs.mpert+1):((j+1)*ffs.mpert)    # mode index i ↔ (m, n) = ((i-1)%mpert + mlow, (i-1)÷mpert + nlow)
+        plasma = ts.toroidal_mode[j+1].plasma
+        plasma.grid_type.index = 24
+        plasma.grid_type.name = "inverse_rhopolnorm_straight_field_line_fourier"
+        plasma.grid.dim1 = sqrt.(pe.psi_grid)
+        plasma.grid.dim2 = ms
+        xi, bn = pe.xi_n_modes[:, cols], pe.b_n_modes[:, cols]
+        plasma.displacement_perpendicular.real = real.(xi)
+        plasma.displacement_perpendicular.imaginary = imag.(xi)
+        plasma.b_field_perturbed.coordinate1.real = real.(bn)
+        plasma.b_field_perturbed.coordinate1.imaginary = imag.(bn)
+    end
+    return nothing
+end
+
+# One toroidal_mode per tearing root, after the ideal entries; tearing is solved for a single n.
+function _write_imas_tearing!(ts, ffs, tearing)
+    (tearing.enabled && !isempty(tearing.gamma) && ffs.npert == 1) || return nothing
+    coupled = tearing.coupled_extraction !== nothing
+    nideal = length(ts.toroidal_mode)
+    resize!(ts.toroidal_mode, nideal + length(tearing.gamma))
+    for (k, (γ, ω)) in enumerate(zip(tearing.gamma, tearing.omega))
+        mode = ts.toroidal_mode[nideal+k]
+        mode.n_tor = ffs.nlow
+        mode.perturbation_type.name = "tearing"
+        mode.perturbation_type.description = coupled ? "coupled tearing root of all rational surfaces" :
+                                             "uncoupled tearing root of one rational surface"
+        coupled || isempty(tearing.rational_q) || (mode.m_pol_dominant = round(tearing.rational_q[k] * ffs.nlow))
+        mode.growthrate = γ
+        mode.frequency = ω / (2π)
+    end
+    return nothing
+end
+
+# ρ_tor,N at the given ψ_N from the dd's active equilibrium slice (the IMAS-fed case), or nothing without one.
+function _rho_tor_norm(dd, psi_norm)
+    isempty(dd.equilibrium.time_slice) && return nothing
+    eqt = dd.equilibrium.time_slice[]
+    p = eqt.profiles_1d
+    (IMASdd.hasdata(p, :psi) && IMASdd.hasdata(p, :rho_tor_norm)) || return nothing
+    g = eqt.global_quantities
+    pn = (p.psi .- g.psi_axis) ./ (g.psi_boundary - g.psi_axis)
+    return linear_interp(pn, p.rho_tor_norm).(psi_norm)
+end
+
+# NTV torque (real part) of each method as a neoclassical core source: density per volume and torque inside ψ_N.
+function _write_imas_torque!(dd, equil, kinetic)
+    for (method, r) in sort!(collect(kinetic.total.method_results); by=first)
+        isempty(r.psi_grid) && continue
+        src = resize!(dd.core_sources.source, :neoclassical, "identifier.description" => "GPEC NTV torque, method $method"; wipe=false)
+        src.identifier.name = "neoclassical"
+        p1 = resize!(src.profiles_1d; wipe=false)
+        rho = _rho_tor_norm(dd, r.psi_grid)
+        rho === nothing || (p1.grid.rho_tor_norm = rho)
+        setproperty!(p1.grid, :psi_norm, copy(r.psi_grid); error_on_missing_coordinates=false)
+        dVdpsi = [equil.profiles.dVdpsi_spline(ψ) for ψ in r.psi_grid]
+        # Without ρ_tor,N (a run not fed from the dd) the profiles are kept on ψ_N alone.
+        setproperty!(p1, :momentum_tor, real.(r.dtdpsi) ./ dVdpsi; error_on_missing_coordinates=false)
+        setproperty!(p1, :torque_tor_inside, real.(r.t_cumulative); error_on_missing_coordinates=false)
+    end
+    return nothing
+end
+
 export main, write_imas
 export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields
 export KineticProfileData, IonSpecies, attach_kinetic_profiles!

test/runtests_imas.jl

@@ -229,4 +229,59 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         @test isempty(dd.mhd_linear.time_slice)
     end
 
+    # Test 8: write_imas — later stages: each n's response block, tearing roots after the ideal entries,
+    # and the torque density, whose volume integral must give the torque inside the edge.
+    @testset "write_imas: response, tearing and torque" begin
+        npsi, mpert, npert = 5, 3, 2
+        psi = collect(range(0.1, 0.9; length=npsi))
+        xi = [complex(ψ, 10j) for ψ in psi, j in 1:(mpert*npert)]
+        dVdpsi(ψ) = 2.0 + ψ                                  # m³ per unit ψ_N
+        dtdpsi = complex.(3.0 .* dVdpsi.(psi))                # uniform torque density of 3 N m⁻²
+        t_cum = complex.([3.0 * (2.0 * (ψ - psi[1]) + (ψ^2 - psi[1]^2) / 2) for ψ in psi])
+        mr = (psi_grid=psi, dtdpsi=dtdpsi, t_cumulative=t_cum, total_torque=t_cum[end])
+        result = (
+            ffs=(integrator=:forward, free_boundary=(et=[0.3 + 0im, 0.5 + 0im, 0.6 + 0im, 0.7 + 0im, 0.8 + 0im, 0.9 + 0im],
+                    n_tor_idx=[0, 1, 0, 1, 0, 1]), numpert_total=mpert * npert, npert=npert, nlow=1, mlow=-1, mhigh=1, mpert=mpert),
+            pe=(psi_grid=psi, xi_n_modes=xi, b_n_modes=2 .* xi),
+            kinetic=(total=(method_results=Dict("fgar" => mr),),),
+            equil=(profiles=(dVdpsi_spline=dVdpsi,),)
+        )
+        dd = IMASdd.dd()
+        dd.global_time = 0.0
+        GeneralizedPerturbedEquilibrium.write_imas(dd, result)
+        tm = dd.mhd_linear.time_slice[1].toroidal_mode
+        @test length(tm) == npert
+        @test tm[2].plasma.displacement_perpendicular.real == real.(xi[:, 4:6])   # the n = 2 block
+        @test tm[1].plasma.b_field_perturbed.coordinate1.imaginary == imag.(2 .* xi[:, 1:3])
+        @test tm[1].plasma.grid.dim1 ≈ sqrt.(psi)
+        @test tm[1].plasma.grid.dim2 == [-1.0, 0.0, 1.0]
+        p1 = dd.core_sources.source[1].profiles_1d[1]
+        @test dd.core_sources.source[1].identifier.index == 401
+        @test dd.core_sources.source[1].identifier.name == "neoclassical"
+        @test p1.momentum_tor ≈ fill(3.0, npsi)
+        @test p1.torque_tor_inside[end] ≈ real(mr.total_torque)
+        @test !IMASdd.hasdata(p1.grid, :rho_tor_norm)        # no equilibrium in the dd
+
+        # With the equilibrium in the dd (COCOS 11), ρ_tor,N is interpolated onto the torque grid
+        dde, _ = make_mock_dd()
+        eqp = dde.equilibrium.time_slice[1].profiles_1d
+        eqp.rho_tor_norm = sqrt.(range(0, 1; length=length(eqp.psi)))
+        GeneralizedPerturbedEquilibrium.write_imas(dde, result)
+        @test dde.core_sources.source[1].profiles_1d[1].grid.rho_tor_norm ≈ sqrt.(psi) rtol = 1e-2
+
+        # Tearing is single n: two uncoupled roots follow the one ideal entry
+        tr = (enabled=true, gamma=[12.0, -3.0], omega=[2π * 100, 0.0], coupled_extraction=nothing, rational_q=[2.0, 3.0])
+        result1 = (ffs=(integrator=:forward, free_boundary=(et=[0.4 + 0im], n_tor_idx=[0]), numpert_total=1, npert=1, nlow=1),
+            tearing=tr)
+        dd1 = IMASdd.dd()
+        GeneralizedPerturbedEquilibrium.write_imas(dd1, result1)
+        tm1 = dd1.mhd_linear.time_slice[1].toroidal_mode
+        @test length(tm1) == 3
+        @test tm1[1].energy_perturbed == 0.4
+        @test [m.growthrate for m in tm1[2:3]] == tr.gamma
+        @test tm1[2].frequency ≈ 100.0
+        @test [m.m_pol_dominant for m in tm1[2:3]] == [2.0, 3.0]
+        @test all(m.perturbation_type.name == "tearing" for m in tm1[2:3])
+    end
+
 end