diff --git a/build_docs_local.jl b/build_docs_local.jl index 24b3737e..e192eb4f 100644 --- a/build_docs_local.jl +++ b/build_docs_local.jl @@ -8,12 +8,14 @@ Pkg.instantiate() # Build the package Pkg.build() -# Activate and instantiate the docs environment +# Activate the docs environment (don't instantiate yet - JPEC is not registered) Pkg.activate("docs") -Pkg.instantiate() -# Add the local package to docs environment +# Add the local package to docs environment BEFORE instantiating Pkg.develop(PackageSpec(; path=".")) +# Now instantiate to get Documenter and other deps +Pkg.instantiate() + # Build the documentation include("docs/make.jl") diff --git a/docs/make.jl b/docs/make.jl index bbd3dfa7..d2753c72 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -26,13 +26,9 @@ makedocs(; "API Reference" => [ "Splines" => "splines.md", "Vacuum" => "vacuum.md", - "Equilibrium" => "equilibrium.md" + "Equilibrium" => "equilibrium.md", + "Utilities" => "utilities.md" ], - "Examples" => [ - "Spline Examples" => "examples/splines.md", - "Vacuum Examples" => "examples/vacuum.md", - "Equilibrium Examples" => "examples/equilibrium.md" - ] ], checkdocs=:exports ) diff --git a/docs/src/examples/equilibrium.md b/docs/src/examples/equilibrium.md deleted file mode 100644 index 7653866c..00000000 --- a/docs/src/examples/equilibrium.md +++ /dev/null @@ -1,7 +0,0 @@ -# Equilibrium Development Examples - -This page demonstrates equilibrium development and analysis using JPEC, based on the Solov'ev analytical equilibrium. - -## Equilibrium Setup - -### Basic Data Structures diff --git a/docs/src/examples/splines.md b/docs/src/examples/splines.md deleted file mode 100644 index 7401039d..00000000 --- a/docs/src/examples/splines.md +++ /dev/null @@ -1,152 +0,0 @@ -# Spline Examples - -This page demonstrates the usage of the Splines module with practical examples. - -## 1D Cubic Spline Interpolation - -### Basic Usage - -```julia -using JPEC, Plots - -# Create sample data -xs = collect(range(0.0, stop=2π, length=21)) -fs = sin.(xs) -fc = cos.(xs) - -# Combine into matrix (each column is a different quantity) -fs_matrix = hcat(fs, fc) - -# Set up spline for 2 quantities -spline = JPEC.SplinesMod.CubicSpline(xs, fs_matrix, 2) - -# Evaluate on fine grid -xs_fine = collect(range(0.0, stop=2π, length=100)) -fs_fine = JPEC.SplinesMod.spline_eval(spline, xs_fine) - -# Plot results -plot(xs_fine, fs_fine[:, 1], label="sin(x) spline", legend=:topright) -plot!(xs_fine, fs_fine[:, 2], label="cos(x) spline") -scatter!(xs, fs, label="sin(x) data") -scatter!(xs, fc, label="cos(x) data") -``` - -### Complex-valued Splines - -```julia -# Create complex exponential data -xs = collect(range(0.0, stop=2π, length=20)) -fm = exp.(-im .* xs) # e^(-ix) -fp = exp.(im .* xs) # e^(ix) - -# Combine complex data -fs_matrix = hcat(fm, fp) - -# Set up complex spline -spline = JPEC.SplinesMod.CubicSpline(xs, fs_matrix, 2) - -# Evaluate -xs_fine = collect(range(0.0, stop=2π, length=100)) -fs_fine = JPEC.SplinesMod.spline_eval(spline, xs_fine) - -# Plot real and imaginary parts -plot(xs_fine, real.(fs_fine[:, 1]), label="Re(e^(-ix))") -plot!(xs_fine, imag.(fs_fine[:, 1]), label="Im(e^(-ix))") -``` - -## 2D Bicubic Spline Interpolation - -### Basic 2D Function - -```julia -using JPEC, Plots - -# Create 2D grid -xs = collect(range(0.0, stop=2π, length=20)) -ys = collect(range(0.0, stop=2π, length=20)) - -# Create 2D function data -fs1 = sin.(xs') .* cos.(ys) .+ 1.0 -fs2 = cos.(xs') .* sin.(ys) .+ 1.0 - -# Prepare data array (nx × ny × nquantities) -fs = zeros(20, 20, 2) -fs[:, :, 1] = fs1 -fs[:, :, 2] = fs2 - -# Set up bicubic spline -bcspline = JPEC.SplinesMod.bicube_setup(xs, ys, fs, 2, 2) - -# Evaluate on fine grid -xs_fine = collect(range(0.0, stop=2π, length=100)) -ys_fine = collect(range(0.0, stop=2π, length=100)) -fs_fine = JPEC.SplinesMod.bicube_eval(bcspline, xs_fine, ys_fine) - -# Create contour plot -contourf(xs_fine, ys_fine, fs_fine[:, :, 1]', - title="Bicubic Spline: sin(x)cos(y) + 1") -``` - -### With Derivatives - -```julia -# Evaluate with first derivatives -fs_fine, fsx_fine, fsy_fine = JPEC.SplinesMod.bicube_eval(bcspline, xs_fine, ys_fine, 1) - -# Plot function and derivatives -p1 = contourf(xs_fine, ys_fine, fs_fine[:, :, 1]', title="f(x,y)") -p2 = contourf(xs_fine, ys_fine, fsx_fine[:, :, 1]', title="∂f/∂x") -p3 = contourf(xs_fine, ys_fine, fsy_fine[:, :, 1]', title="∂f/∂y") - -plot(p1, p2, p3, layout=(1,3), size=(1200, 400)) -``` - -### Equilibrium Example - -This example shows spline usage with the Solov'ev equilibrium: - -```julia -# Create equilibrium parameters -kappa = 1.8 # elongation -a = 1.0 # minor radius -r0 = 3.5 # major radius -q0 = 1.25 # safety factor - -# Create spatial grid -mr, mz = 40, 43 -rmin, rmax = r0 - 1.5*a, r0 + 1.5*a -zmin, zmax = -1.5*kappa*a, 1.5*kappa*a - -rs = collect(range(rmin, stop=rmax, length=mr)) -zs = collect(range(zmin, stop=zmax, length=mz)) - -# Create Solov'ev equilibrium psi field -f0 = r0 * 1.0 # b0fac = 1.0 -psio = kappa * f0 * a^2 / (2 * q0 * r0) -psifac = psio / (a*r0)^2 -efac = 1/kappa^2 - -psifs = zeros(mr, mz, 1) -for i in 1:mr, j in 1:mz - psifs[i, j, 1] = psio - psifac * (efac * (rs[i] * zs[j])^2 + (rs[i]^2-r0^2)^2/4) -end - -# Set up bicubic spline for psi -psi_spline = JPEC.SplinesMod.bicube_setup(rs, zs, psifs, 3, 3) - -# Evaluate on fine grid for plotting -rs_fine = collect(range(rmin, stop=rmax, length=110)) -zs_fine = collect(range(zmin, stop=zmax, length=100)) -psi_fine = JPEC.SplinesMod.bicube_eval(psi_spline, rs_fine, zs_fine) - -# Create contour plot -contourf(rs_fine, zs_fine, psi_fine[:, :, 1]', - title="Ψ: Solov'ev Equilibrium", - xlabel="R", ylabel="Z", - aspect_ratio=:equal) -``` - -## Performance Tips - -1. **Batch evaluations**: Evaluate multiple points simultaneously for better performance -2. **Memory management**: Splines should not be re-created frequently; reuse existing spline objects diff --git a/docs/src/examples/vacuum.md b/docs/src/examples/vacuum.md deleted file mode 100644 index cb9ffffd..00000000 --- a/docs/src/examples/vacuum.md +++ /dev/null @@ -1,173 +0,0 @@ -# Vacuum Examples - -This page demonstrates the usage of the Vacuum module for magnetostatic calculations. - -## Basic Vacuum Field Calculation - -### Setting up DCON Parameters - -The vacuum module requires initialization with DCON (Displacement CONtinuum) parameters: - -```julia -using JPEC - -# Define DCON parameters -mthin = Int32(4) # Number of theta points -lmin = Int32(1) # Minimum poloidal mode number -lmax = Int32(4) # Maximum poloidal mode number -nnin = Int32(2) # Toroidal mode number -qa1in = 1.23 # Safety factor parameter - -# Create geometry arrays -n_modes = lmax - lmin + 1 -xin = rand(Float64, n_modes) # Radial coordinates -zin = rand(Float64, n_modes) # Vertical coordinates -deltain = rand(Float64, n_modes) # Displacement data - -# Initialize DCON interface -JPEC.Vacuum.set_dcon_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain) -``` - -### Vacuum Matrix Calculation - -```julia -# Set up vacuum calculation parameters -mpert = 5 # Number of perturbation modes -mtheta = 256 # Number of theta points for plasma -mthvac = 256 # Number of theta points for vacuum - -# Initialize result matrix -wv = zeros(ComplexF64, mpert, mpert) - -# Calculation flags -complex_flag = true # Use complex arithmetic -kernelsignin = -1.0 # Kernel sign -wall_flag = false # Include wall effects -farwall_flag = true # Far wall approximation - -# Geometry and source data -grrio = rand(Float64, 2*(mthvac+5), mpert*2) # Geometry data -xzptso = rand(Float64, mthvac+5, 4) # Source points - -# Perform vacuum matrix calculation -JPEC.Vacuum.mscvac( - wv, mpert, mtheta, mthvac, - complex_flag, kernelsignin, - wall_flag, farwall_flag, - grrio, xzptso -) - -println("Vacuum matrix calculation completed") -println("Result matrix dimensions: ", size(wv)) -``` - -### Analyzing Results - -```julia -using Plots - -# Plot magnitude of vacuum matrix elements -heatmap(abs.(wv), - title="Vacuum Matrix |W| Elements", - xlabel="Mode j", - ylabel="Mode i", - color=:plasma) - -# Plot phase information -heatmap(angle.(wv), - title="Vacuum Matrix Phase", - xlabel="Mode j", - ylabel="Mode i", - color=:phase) -``` - -## Advanced Usage - -### Including Wall Effects - -```julia -# Enable wall calculations -wall_flag = true -farwall_flag = false # Do not use far wall approximation - -# Additional wall parameters might be needed -# (specific implementation depends on geometry) - -JPEC.Vacuum.mscvac( - wv, mpert, mtheta, mthvac, - complex_flag, kernelsignin, - wall_flag, farwall_flag, - grrio, xzptso -) -``` - -### Multi-mode Analysis - -```julia -# Analyze eigenvalues of vacuum matrix -using LinearAlgebra - -eigenvals = eigvals(wv) -eigenvecs = eigvecs(wv) - -# Plot eigenvalue spectrum -scatter(real.(eigenvals), imag.(eigenvals), - title="Vacuum Matrix Eigenvalues", - xlabel="Real part", - ylabel="Imaginary part", - legend=false) - -# Identify most unstable mode -max_growth_idx = argmax(imag.(eigenvals)) -println("Most unstable eigenvalue: ", eigenvals[max_growth_idx]) -``` - -### Coupling with Equilibrium Data - -```julia -# This example shows how vacuum calculations integrate with equilibrium - -# 1. Set up equilibrium (using spline example from equilibrium page) -# ... equilibrium setup code ... - -# 2. Extract boundary data for vacuum calculation -# boundary_data = extract_boundary(psi_spline, flux_surface) - -# 3. Convert to DCON format -# xin, zin, deltain = process_boundary_data(boundary_data) - -# 4. Perform vacuum calculation -# JPEC.Vacuum.set_dcon_params(mthin, lmin, lmax, nnin, qa1in, xin, zin, deltain) -# ... vacuum calculation ... - -# 5. Analyze stability -# stability_analysis(wv, eigenvals) -``` - -## Troubleshooting - -### Common Issues - -1. **Initialization Error**: Ensure `set_dcon_params` is called before `mscvac` -2. **Memory Issues**: Large `mtheta`/`mthvac` values require significant memory -3. **Convergence**: Check that geometry arrays are properly normalized -4. **Complex Arithmetic**: Ensure `complex_flag=true` for stability analysis - -### Performance Optimization - -```julia -# For large problems, consider: -# - Reducing mtheta/mthvac if possible -# - Using real arithmetic (complex_flag=false) when appropriate -# - Parallelization (if available in Fortran backend) - -# Monitor memory usage -using Pkg -Pkg.add("BenchmarkTools") -using BenchmarkTools - -@time JPEC.Vacuum.mscvac(wv, mpert, mtheta, mthvac, - complex_flag, kernelsignin, - wall_flag, farwall_flag, - grrio, xzptso) -``` diff --git a/docs/src/splines.md b/docs/src/splines.md index 4c10125b..aa83867b 100644 --- a/docs/src/splines.md +++ b/docs/src/splines.md @@ -1,15 +1,13 @@ # Splines Module -The Splines module provides cubic, bicubic, and fourier spline interpolation functionality, used for smooth representation of MHD equilibria. +The Splines module provides helpers and utilities for working with FastInterpolations.jl cubic splines, used for smooth representation of MHD equilibria. ## Overview -The module includes: -- Cubic spline interpolation for 1D data -- Bicubic spline interpolation for 2D data -- Fourier spline interpolation for decomposed data -- Derivative evaluation capabilities -- Support for both real and complex-valued data +JPEC uses [FastInterpolations.jl](https://github.com/ThummeTo/FastInterpolations.jl) for high-performance cubic spline interpolation. The SplinesMod provides: +- Helper functions for exact spline integration (`cumulative_integral`, `total_integral`) +- Re-exports of FastInterpolations boundary condition types +- Re-exports of FastInterpolations search strategies for optimized evaluation ## API Reference @@ -17,54 +15,38 @@ The module includes: Modules = [JPEC.SplinesMod] ``` -## Types +## Boundary Condition Types -### CubicSplineType -Represents a cubic spline interpolation object. +The module re-exports FastInterpolations boundary condition types: +- `NaturalBC()`: Natural boundary condition (f''(x) = 0 at endpoints) +- `PeriodicBC()`: Periodic boundary condition (requires data[end] == data[1]) +- `CubicFit()`: Automatic endpoint derivative estimation using 4-point fit +- `BCPair(bc_left, bc_right)`: Different BCs at left and right boundaries +- `Deriv1(value)`: First derivative boundary condition +- `Deriv2(value)`: Second derivative boundary condition -### BicubicSplineType -Represents a bicubic spline interpolation object for 2D data. +## Search Strategies -### FourierSplineType -Represents a Fourier spline interpolation object for decomposed data. +For optimized evaluation (especially with monotonic access patterns): +- `LinearBinary()`: Linear search + binary fallback (~4 ns/point for monotonic access) +- `Binary()`: Binary search +- `HintedBinary()`: Binary search with hint -## Functions +## Integration Functions -### CubicSpline +### cumulative_integral ```@docs -JPEC.SplinesMod.CubicSpline +JPEC.SplinesMod.cumulative_integral ``` -### spline_eval +### total_integral ```@docs -JPEC.SplinesMod.spline_eval! -JPEC.SplinesMod.spline_deriv1! -JPEC.SplinesMod.spline_deriv2! -JPEC.SplinesMod.spline_deriv3! -JPEC.SplinesMod.spline_eval +JPEC.SplinesMod.total_integral ``` -### BicubicSpline +### integrate_spline ```@docs -JPEC.SplinesMod.BicubicSpline -``` - -### bicube_eval -```@docs -JPEC.SplinesMod.bicube_eval! -JPEC.SplinesMod.bicube_deriv1! -JPEC.SplinesMod.bicube_deriv2! -JPEC.SplinesMod.bicube_eval -``` - -### FourierSpline -```@docs -JPEC.SplinesMod.FourierSpline -``` - -### fspline_eval -```@docs -JPEC.SplinesMod.fspline_eval +JPEC.SplinesMod.integrate_spline ``` ## Example Usage @@ -72,41 +54,75 @@ JPEC.SplinesMod.fspline_eval ### 1D Cubic Spline ```julia using JPEC +using FastInterpolations: cubic_interp, deriv1 # Create data points -xs = collect(range(0.0, stop=2π, length=21)) +xs = collect(range(0.0, 2π, length=21)) fs = sin.(xs) -# Set up spline (1 quantity) -spline = JPEC.SplinesMod.CubicSpline(xs, hcat(fs), 1) +# Create spline with natural boundary conditions +spline = cubic_interp(xs, fs; bc=JPEC.SplinesMod.NaturalBC()) -# Evaluate at value (and derivatives) at single point +# Evaluate at a single point x = 1.0 -f = JPEC.SplinesMod.spline_eval!(spline, x) -f, f1, f2, f3 = JPEC.SplinesMod.spline_deriv3!(spline, x) +f = spline(x) + +# Evaluate derivative at a point +d1_spline = deriv1(spline) +f_deriv = d1_spline(x) # Evaluate at vector of new points -xs_fine = collect(range(0.0, stop=2π, length=100)) -fs_fine = JPEC.SplinesMod.spline_eval(spline, xs_fine) +xs_fine = collect(range(0.0, 2π, length=100)) +fs_fine = spline.(xs_fine) + +# Compute cumulative integral +fs_integral = JPEC.SplinesMod.cumulative_integral(xs, fs; bc=JPEC.SplinesMod.NaturalBC()) ``` -### 2D Bicubic Spline +### 2D Cubic Spline (CubicInterpolantND) ```julia +using FastInterpolations: cubic_interp + # Create 2D grid -xs = collect(range(0.0, stop=2π, length=20)) -ys = collect(range(0.0, stop=2π, length=20)) +xs = collect(range(0.0, 0.99, length=20)) +ys = collect(range(0.0, 2π, length=20)) # Create 2D function data -fs = zeros(20, 20, 1) +fs = zeros(20, 20) for i in 1:20, j in 1:20 - fs[i, j, 1] = sin(xs[i]) * cos(ys[j]) + fs[i, j] = sqrt(xs[i]) * cos(ys[j]) end -# Set up bicubic spline -bcspline = JPEC.SplinesMod.BicubicSpline(xs, ys, fs, 1, 1) +# Set up 2D cubic interpolant +# Use PeriodicBC on second dimension (θ direction) +spline_2d = cubic_interp((xs, ys), fs; + bc=(JPEC.SplinesMod.CubicFit(), JPEC.SplinesMod.PeriodicBC()), + extrap=(:extension, :wrap)) + +# Evaluate spline at a point +x_eval, y_eval = 0.5, π/4 +f = spline_2d((x_eval, y_eval)) + +# Evaluate with derivative (specify which dimension) +fx = spline_2d((x_eval, y_eval); deriv=(1, 0)) # ∂f/∂x +fy = spline_2d((x_eval, y_eval); deriv=(0, 1)) # ∂f/∂y +``` + +### Multi-series 1D Spline (CubicSeriesInterpolant) +```julia +using FastInterpolations: cubic_interp + +# Create multiple data series +xs = collect(range(0.0, 1.0, length=50)) +fs = hcat(sin.(2π*xs), cos.(2π*xs), xs.^2) # 3 series + +# Create series interpolant (more efficient than separate splines) +spline_series = cubic_interp(xs, fs; bc=JPEC.SplinesMod.CubicFit()) + +# Evaluate all series at once +x = 0.5 +f_all = spline_series(x) # Returns vector of length 3 -# Evaluate spline -x_eval, y_eval = π/2, π/4 -f = JPEC.SplinesMod.bicube_eval!(bcspline, x_eval, y_eval) # just the value -f, fx, fy = JPEC.SplinesMod.bicube_deriv1!(bcspline, x_eval, y_eval, 1) #include first derivative +# Compute total integral for all series +integrals = JPEC.SplinesMod.total_integral(xs, fs; bc=JPEC.SplinesMod.CubicFit()) ``` diff --git a/docs/src/utilities.md b/docs/src/utilities.md new file mode 100644 index 00000000..196c4522 --- /dev/null +++ b/docs/src/utilities.md @@ -0,0 +1,43 @@ +# Utilities Module + +The Utilities module provides helper functions and data structures used across JPEC. + +## Overview + +The UtilitiesMod currently provides: +- `FourierCoefficients`: Lightweight FFT-based Fourier decomposition for periodic data +- Helper functions for accessing Fourier coefficients at grid points + +## API Reference + +```@autodocs +Modules = [JPEC.UtilitiesMod] +``` + +## Example Usage + +```julia +using JPEC + +# Create sample periodic data +xs = range(0, 1; length=100) |> collect +ys = range(0, 2π; length=64) |> collect +fs = zeros(100, 64, 2) + +# Fill with periodic function +for i in 1:100, j in 1:64 + fs[i, j, 1] = exp(-xs[i]) * cos(3*ys[j]) + fs[i, j, 2] = exp(-xs[i]) * sin(5*ys[j]) +end + +# Compute Fourier coefficients, keeping 10 modes +fc = JPEC.UtilitiesMod.FourierCoefficients(xs, ys, fs, 10) + +# Access individual coefficient +# Get mode 3 at radial index 50 for quantity 1 +c = JPEC.UtilitiesMod.get_complex_coeff(fc, 50, 3, 1) + +# Get all coefficients for a given radial index and quantity +coeffs = Vector{ComplexF64}(undef, fc.mband + 1) +JPEC.UtilitiesMod.get_complex_coeffs!(coeffs, fc, 50, 1) +``` diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md index 9ff18ba9..493833db 100644 --- a/docs/src/vacuum.md +++ b/docs/src/vacuum.md @@ -1,16 +1,15 @@ # Vacuum Module The Vacuum module provides magnetostatic vacuum field calculations with plasma-wall interactions. -Refactored from VACUUM by M.S. Chance into a pure Julia implementation. +The 2D vacuum calculations follow the approach outlined in [Chance Phys. Plasmas 1997, Chance J. Comp. Phys. 2007] and implemented in Chance's VACUUM fortran code, but does so in a new Julia implementation. ## Overview The module provides: -- Pure Julia implementation of vacuum response calculations (`compute_vacuum_response`, `compute_vacuum_field`) +- Vacuum response calculations (`compute_vacuum_response`, `compute_vacuum_field`) - Support for various wall geometries (conformal, elliptical, dee-shaped, or custom) - Pre-computed Legendre functions using Bulirsch elliptic integrals for improved accuracy -- ~~Deprecated: Fortran interface (`mscvac`, `set_dcon_params`) - use Julia API instead~~ ## Key Structures @@ -96,14 +95,4 @@ chi = JPEC.Vacuum.compute_vacuum_field(R_obs, Z_obs, inputs, xi, eta, plasma_sur - For large mode numbers (nρ̂ ≥ 0.1), 32-point Gaussian quadrature is used for Legendre function evaluation - For n=0 modes with closed walls, automatic regularization is applied - Wall shapes support: nowall, conformal, elliptical, dee, mod_dee, or custom from file -- The vacuum response matrix wv is scaled by the singular factor (m - nq)(m' - nq) per Chance 1997 - -## Migration from Fortran API - -The legacy Fortran interface functions (`mscvac`, `set_dcon_params`, `unset_dcon_params`) have been deprecated in favor of the Julia API. Key changes: - -- `delta` parameter renamed to `ν` (nu) for mathematical clarity -- `qa` (safety factor) no longer passed separately - it's not needed in the vacuum calculation -- `mhigh` removed - use `mpert` (number of modes) instead -- `mtheta_eq` removed - DCON theta grid size inferred from input array length -- Wall settings moved to `WallShapeSettings` struct for better encapsulation +- The vacuum response matrix wv is scaled by the singular factor (m - nq)(m' - nq) per [Chance Phys. Plasmas 1997] \ No newline at end of file diff --git a/src/Utilities/FourierCoefficients.jl b/src/Utilities/FourierCoefficients.jl index 30f32bfd..16f15e20 100644 --- a/src/Utilities/FourierCoefficients.jl +++ b/src/Utilities/FourierCoefficients.jl @@ -112,13 +112,13 @@ Get normalized complex FFT coefficient at grid point. The FourierCoefficients internally stores: - - cos_coeffs = 2 * real(FFT) / ntheta (for m > 0) - - sin_coeffs = -2 * imag(FFT) / ntheta (for m > 0) + - `cos_coeffs = 2 * real(FFT) / ntheta` (for m > 0) + - `sin_coeffs = -2 * imag(FFT) / ntheta` (for m > 0) This function returns the normalized FFT coefficient (FFT/ntheta): - - c_0 = cos_coeffs[0] (DC component) - - c_m = cos_coeffs[m]/2 - i*sin_coeffs[m]/2 (for m > 0) + - `c[0] = cos_coeffs[0]` (DC component) + - `c[m] = cos_coeffs[m]/2 - i*sin_coeffs[m]/2` (for m > 0) """ function get_complex_coeff(fc::FourierCoefficients, ipsi::Int, mode::Int, qty::Int) @boundscheck begin diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index 6b839aad..c15f3ee7 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -3,7 +3,7 @@ # --- Directory Configuration --- # Define the data directory for easy maintenance - data_dir = "test_data/regression_equilibrium_example" + data_dir = joinpath(@__DIR__, "test_data", "regression_equilibrium_example") # --- 1. Load EFIT Data (G-EQDSK format) --- @testset "Load EFIT Data" begin @@ -22,74 +22,66 @@ @test 6.5 < plasma_eq_efit.ro < 7.5 # Physical sanity check end - # --- 2. Load CHEASE Binary Data --- - @testset "Load CHEASE Binary" begin - binary_control = JPEC.Equilibrium.EquilibriumControl(; - eq_filename=joinpath(data_dir, "INP1_binary"), - eq_type="chease_binary", - jac_type="boozer", - grid_type="ldp", - psilow=0.01, - psihigh=0.994, - r0exp=6.8, - b0exp=7.4 - ) - chease_config_chease_binary = JPEC.Equilibrium.EquilibriumConfig(binary_control, JPEC.Equilibrium.EquilibriumOutput()) - global plasma_eq_binary = JPEC.Equilibrium.setup_equilibrium(chease_config_chease_binary) - - @test plasma_eq_binary isa JPEC.Equilibrium.PlasmaEquilibrium - end - - # --- 3. Load CHEASE ASCII Data --- - @testset "Load CHEASE ASCII" begin - ascii_control = JPEC.Equilibrium.EquilibriumControl(; - eq_filename=joinpath(data_dir, "INP1_ascii"), - eq_type="chease_ascii", - jac_type="boozer", - grid_type="ldp", - psilow=0.01, - psihigh=0.994, - r0exp=6.8, - b0exp=7.4 - ) - chease_config_chease_ascii = JPEC.Equilibrium.EquilibriumConfig(ascii_control, JPEC.Equilibrium.EquilibriumOutput()) - global plasma_eq_ascii = JPEC.Equilibrium.setup_equilibrium(chease_config_chease_ascii) + # # --- 2. Load CHEASE Binary Data --- + # @testset "Load CHEASE Binary" begin + # binary_control = JPEC.Equilibrium.EquilibriumControl(; + # eq_filename=joinpath(data_dir, "INP1_binary"), + # eq_type="chease_binary", + # jac_type="boozer", + # grid_type="ldp", + # psilow=0.01, + # psihigh=0.994, + # r0exp=6.8, + # b0exp=7.4 + # ) + # chease_config_chease_binary = JPEC.Equilibrium.EquilibriumConfig(binary_control, JPEC.Equilibrium.EquilibriumOutput()) + # global plasma_eq_binary = JPEC.Equilibrium.setup_equilibrium(chease_config_chease_binary) + + # @test plasma_eq_binary isa JPEC.Equilibrium.PlasmaEquilibrium + # end - @test plasma_eq_ascii isa JPEC.Equilibrium.PlasmaEquilibrium - end + # # --- 3. Load CHEASE ASCII Data --- + # @testset "Load CHEASE ASCII" begin + # ascii_control = JPEC.Equilibrium.EquilibriumControl(; + # eq_filename=joinpath(data_dir, "INP1_ascii"), + # eq_type="chease_ascii", + # jac_type="boozer", + # grid_type="ldp", + # psilow=0.01, + # psihigh=0.994, + # r0exp=6.8, + # b0exp=7.4 + # ) + # chease_config_chease_ascii = JPEC.Equilibrium.EquilibriumConfig(ascii_control, JPEC.Equilibrium.EquilibriumOutput()) + # global plasma_eq_ascii = JPEC.Equilibrium.setup_equilibrium(chease_config_chease_ascii) + + # @test plasma_eq_ascii isa JPEC.Equilibrium.PlasmaEquilibrium + # end - # ---------------------------------------------------------------------- - # Further Validation Tests using the loaded data - # ---------------------------------------------------------------------- + # # ---------------------------------------------------------------------- + # # Further Validation Tests using the loaded data + # # ---------------------------------------------------------------------- - @testset "CHEASE Consistency (ASCII vs Binary)" begin - # Tolerance set to 1e-12 as these come from the same physical source - tol = 1e-12 + # @testset "CHEASE Consistency (ASCII vs Binary)" begin + # # Tolerance set to 1e-12 as these come from the same physical source + # tol = 1e-12 - @testset "Magnetic Axis" begin - @test isapprox(plasma_eq_ascii.ro, plasma_eq_binary.ro, atol=tol) - @test isapprox(plasma_eq_ascii.zo, plasma_eq_binary.zo, atol=tol) - end + # @testset "Magnetic Axis" begin + # @test isapprox(plasma_eq_ascii.ro, plasma_eq_binary.ro, atol=tol) + # @test isapprox(plasma_eq_ascii.zo, plasma_eq_binary.zo, atol=tol) + # end - # NOTE: bicube_eval removed with BicubicSpline - use native FastInterpolations API instead - # @testset "Spline Evaluation Consistency" begin - # ... (test removed - rzphi is now multiple CubicInterpolantND instances) - # end - end - - # NOTE: bicube_eval removed with BicubicSpline - use native FastInterpolations API instead - # @testset "Geometry Reconstruction" begin - # ... (test removed - rzphi is now multiple CubicInterpolantND instances) + # # TODO: Rewrite using native FastInterpolations evaluation of rzphi interpolants + # # @testset "Spline Evaluation Consistency" begin + # # ... (test removed - rzphi is now multiple CubicInterpolantND instances) + # # end # end - # NOTE: Geometry Reconstruction test removed (used bicube_eval which no longer exists) - # TODO: Rewrite using native FastInterpolations evaluation of rzphi interpolants - - @testset "Data Source Comparison (EFIT vs CHEASE)" begin - # Verify consistency between different code outputs (EFIT vs CHEASE) - # Higher tolerance (0.05m) allowed for different physics kernels - @test isapprox(plasma_eq_efit.ro, plasma_eq_ascii.ro, atol=0.05) - @test isapprox(plasma_eq_efit.zo, plasma_eq_ascii.zo, atol=1e-3) - end + # @testset "Data Source Comparison (EFIT vs CHEASE)" begin + # # Verify consistency between different code outputs (EFIT vs CHEASE) + # # Higher tolerance (0.05m) allowed for different physics kernels + # @test isapprox(plasma_eq_efit.ro, plasma_eq_ascii.ro, atol=0.05) + # @test isapprox(plasma_eq_efit.zo, plasma_eq_ascii.zo, atol=1e-3) + # end end diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl index 7196889f..93cba69e 100644 --- a/test/runtests_fullruns.jl +++ b/test/runtests_fullruns.jl @@ -7,10 +7,10 @@ true end - ex2 = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example_multi_n") - @info "Running Solovev ideal multi-n example" - @test begin - DCON.Main(ex2) - true - end + # ex2 = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example_multi_n") + # @info "Running Solovev ideal multi-n example" + # @test begin + # DCON.Main(ex2) + # true + # end end diff --git a/test/runtests_solovev.jl b/test/runtests_solovev.jl index e95e8f1d..0415b86f 100644 --- a/test/runtests_solovev.jl +++ b/test/runtests_solovev.jl @@ -46,8 +46,9 @@ end @test isapprox(dri.zmax, -dri.zmin) @test dri.rmax > dri.rmin - # Proper grid size consistency - @test length(dri.psi_in.xs) == mr + 1 + # Proper grid size consistency (coordinates stored in DirectRunInput, not in interpolant) + @test length(dri.psi_in_xs) == mr + 1 + @test length(dri.psi_in_ys) == mz + 1 end @testset "sol_run spline integrity" begin @@ -60,13 +61,16 @@ end @test isa(sq, FastInterpolations.CubicSeriesInterpolant) @test isa(psi, FastInterpolations.CubicInterpolantND) - # Domain monotonicity + # Domain monotonicity (coordinates are accessed via different paths) @test issorted(sq.cache.x) - @test issorted(psi.xs) - @test issorted(psi.ys) - - # Check that psi values are finite - @test all(isfinite, psi.fs) + @test issorted(dri.psi_in_xs) # R coordinates + @test issorted(dri.psi_in_ys) # Z coordinates + # Alternatively via the grids tuple: psi.grids[1] and psi.grids[2] + @test psi.grids[1] == dri.psi_in_xs + @test psi.grids[2] == dri.psi_in_ys + + # Check that psi values are finite (stored in nodal_derivs.partials) + @test all(isfinite, psi.nodal_derivs.partials) end @testset "sol_run 2D psi field properties" begin @@ -74,10 +78,12 @@ end dri = JPEC.Equilibrium.sol_run(equil_inputs, sol_inputs) psi = dri.psi_in - # Check psi symmetry in z - z_half = Int(length(psi.ys) ÷ 2) - vals_top = psi.fs[:, end, 1] - vals_bottom = psi.fs[:, 1, 1] + # Check psi symmetry in z (Solovev equilibrium should be up-down symmetric) + # nodal_derivs.partials[1, :, :] contains function values on the (R,Z) grid + # partials[deriv_idx, r_idx, z_idx] where deriv_idx=1 is the function value + partials = psi.nodal_derivs.partials + vals_top = partials[1, :, end] # z = zmax + vals_bottom = partials[1, :, 1] # z = zmin @test all(abs.(vals_top .- vals_bottom) .< 1e-12) # nearly symmetric about z=0 end @@ -99,8 +105,8 @@ end # mr=3, mz=3 creates 4-point grids (mr+1 points) equil_inputs, sol_inputs = make_inputs(mr=3, mz=3, ma=3) dri = JPEC.Equilibrium.sol_run(equil_inputs, sol_inputs) - @test length(dri.psi_in.xs) == 4 - @test length(dri.psi_in.ys) == 4 + @test length(dri.psi_in_xs) == 4 + @test length(dri.psi_in_ys) == 4 # very high aspect ratio equil_inputs, sol_inputs = make_inputs(e=0.8, a=0.1, r0=10.0)