From d342c41d622eaa64f816d3f5f5a7feb4a412e0d2 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Tue, 6 Jan 2026 13:07:11 -0500 Subject: [PATCH 1/2] VACUUM - IMPROVEMENT - major updates to structure of kernel function and commenting lines for clarity --- src/Vacuum/Vacuum.jl | 123 +++---- src/Vacuum/VacuumInternals.jl | 645 ++++++++++++++++------------------ 2 files changed, 360 insertions(+), 408 deletions(-) diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 081821df2..e4d9da32e 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -24,14 +24,14 @@ Initialize DCON parameters for vacuum field calculations. # Arguments -- `mtheta`: Number of theta grid points (Integer) -- `lmin`: Minimum poloidal mode number (Integer) -- `lmax`: Maximum poloidal mode number (Integer) -- `nnin`: Toroidal mode number (Integer) -- `qa1in`: Safety factor parameter (Float64) -- `xin`: Vector of radial coordinates at plasma boundary (Vector{Float64}) -- `zin`: Vector of vertical coordinates at plasma boundary (Vector{Float64}) -- `deltain`: Vector of displacement values (Vector{Float64}) + - `mtheta`: Number of theta grid points (Integer) + - `lmin`: Minimum poloidal mode number (Integer) + - `lmax`: Maximum poloidal mode number (Integer) + - `nnin`: Toroidal mode number (Integer) + - `qa1in`: Safety factor parameter (Float64) + - `xin`: Vector of radial coordinates at plasma boundary (Vector{Float64}) + - `zin`: Vector of vertical coordinates at plasma boundary (Vector{Float64}) + - `deltain`: Vector of displacement values (Vector{Float64}) # Note @@ -75,8 +75,8 @@ and resets the internal DCON state for future vacuum calculations. # Notes -- Must be called after `set_dcon_params` if you want to reset the DCON memory. -- No arguments are required. + - Must be called after `set_dcon_params` if you want to reset the DCON memory. + - No arguments are required. # Example @@ -99,22 +99,22 @@ Compute the vacuum response matrix for magnetostatic perturbations. # Arguments -- `wv`: Pre-allocated complex matrix (mpert × mpert) to store vacuum response (Array{ComplexF64,2}) -- `mpert`: Number of perturbation modes (Integer) -- `mtheta`: Number of theta grid points for plasma (Integer) -- `mtheta_vacuum`: Number of theta grid points for vacuum region (Integer) -- `complex_flag`: Whether to use complex arithmetic (Bool) -- `kernelsignin`: Sign convention for vacuum kernels (Float64, typically -1.0) -- `wall_flag`: Whether to include an externally defined wall shape (Bool) -- `farwall_flag`: Whether to use far-wall approximation (Bool) -- `grrio`: Green's function data (Array{Float64,2}) -- `xzptso`: Source point coordinates (Array{Float64,2}) -- `op_ahgfile`: Optional communication file for when set_dcon_params is not called (String or Nothing) + - `wv`: Pre-allocated complex matrix (mpert × mpert) to store vacuum response (Array{ComplexF64,2}) + - `mpert`: Number of perturbation modes (Integer) + - `mtheta`: Number of theta grid points for plasma (Integer) + - `mtheta_vacuum`: Number of theta grid points for vacuum region (Integer) + - `complex_flag`: Whether to use complex arithmetic (Bool) + - `kernelsignin`: Sign convention for vacuum kernels (Float64, typically -1.0) + - `wall_flag`: Whether to include an externally defined wall shape (Bool) + - `farwall_flag`: Whether to use far-wall approximation (Bool) + - `grrio`: Green's function data (Array{Float64,2}) + - `xzptso`: Source point coordinates (Array{Float64,2}) + - `op_ahgfile`: Optional communication file for when set_dcon_params is not called (String or Nothing) # Returns -- Modifies `wv` in-place with the computed vacuum response matrix -- Returns the modified `wv` matrix + - Modifies `wv` in-place with the computed vacuum response matrix + - Returns the modified `wv` matrix # Note @@ -201,28 +201,28 @@ end """ compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) -Compute the vacuum response matrix using provided vacuum inputs. +Compute the vacuum response matrix using provided vacuum inputs. -This is the pure Julia implementation that replaces the Fortran `mscvac` function. +This is the pure Julia implementation that replaces the Fortran `mscvac` function. It returns the relevant arrays: `wv`, `grri`, and `xzpts`. # Arguments -- `inputs::VacuumInput`: Struct containing vacuum calculation parameters including mode numbers, - grid resolution, toroidal mode number, and plasma boundary information. -- `wall_settings::WallShapeSettings`: Struct specifying the wall geometry configuration. + - `inputs::VacuumInput`: Struct containing vacuum calculation parameters including mode numbers, + grid resolution, toroidal mode number, and plasma boundary information. + - `wall_settings::WallShapeSettings`: Struct specifying the wall geometry configuration. # Returns -- `wv`: Complex vacuum response matrix (mpert × mpert) relating plasma perturbations to vacuum response -- `grri`: Green's function response matrix (2*mtheta × 2*mpert) in Fourier space -- `xzpts`: Coordinate array (mtheta × 4) containing [R_plasma, Z_plasma, R_wall, Z_wall] + - `wv`: Complex vacuum response matrix (mpert × mpert) relating plasma perturbations to vacuum response + - `grri`: Green's function response matrix (2*mtheta × 2*mpert) in Fourier space + - `xzpts`: Coordinate array (mtheta × 4) containing [R_plasma, Z_plasma, R_wall, Z_wall] # Notes -- This function initializes the plasma surface and wall geometry internally -- The vacuum response includes plasma-plasma and plasma-wall coupling effects -- For n=0 modes with closed walls, a regularization factor is added to prevent singularities + - This function initializes the plasma surface and wall geometry internally + - The vacuum response includes plasma-plasma and plasma-wall coupling effects + - For n=0 modes with closed walls, a regularization factor is added to prevent singularities """ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) @@ -236,7 +236,7 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe # Plasma–Plasma block j1, j2 = 1, 1 - kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, plasma_surf.x, plasma_surf.z, j1, j2, 1, n) + kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, plasma_surf.x, plasma_surf.z, j1, j2, n) # Fourier transform plasma-plasma block fourier_transform!(grri, greenfunction_temp, plasma_surf.cslth, 0, 0) @@ -245,15 +245,15 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe !wall.nowall && begin # Plasma–Wall block j1, j2 = 1, 2 - kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, wall.x, wall.z, j1, j2, 0, n) + kernel!(grad_greenfunction_mat, greenfunction_temp, plasma_surf.x, plasma_surf.z, wall.x, wall.z, j1, j2, n) # Wall–Wall block j1, j2 = 2, 2 - kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, wall.x, wall.z, j1, j2, 0, n) + kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, wall.x, wall.z, j1, j2, n) # Wall–Plasma block j1, j2 = 2, 1 - kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, plasma_surf.x, plasma_surf.z, j1, j2, 1, n) + kernel!(grad_greenfunction_mat, greenfunction_temp, wall.x, wall.z, plasma_surf.x, plasma_surf.z, j1, j2, n) # Fourier transform wall blocks into grri fourier_transform!(grri, greenfunction_temp, plasma_surf.cslth, mtheta, 0) @@ -274,7 +274,7 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe (kernelsign < 0) && begin grad_greenfunction_mat .*= kernelsign # Account for factor of 2 in diagonal terms in eq. 90 of Chance - for i in 1:2 * mtheta + for i in 1:(2*mtheta) grad_greenfunction_mat[i, i] += 2.0 end end @@ -313,6 +313,8 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe end wv = complex.(vacmat, vacmti) + display(wv) + # Create xzpts array xzpts = zeros(Float64, inputs.mtheta, 4) @views xzpts[:, 1] .= plasma_surf.x @@ -323,40 +325,41 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe end """ - compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall::WallGeometry, + compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall::WallGeometry, Bn::Vector{<:Number}, R_grid::AbstractVector, Z_grid::AbstractVector) Calculate the perturbed magnetic field in the vacuum region resulting from a normal magnetic field perturbation (`Bn`) at the plasma surface. Replaces `mscfld` from Fortran. This function orchestrates the vacuum field calculation by: -1. Calling `vaccal!` to compute the vacuum response kernel (`grri`) -2. Defining a grid of points (`R_grid`, `Z_grid`) where the field is to be calculated -3. Calling `_pickup_field` to compute the magnetic field components on that grid using the kernel - and the source perturbation `Bn` + + 1. Calling `vaccal!` to compute the vacuum response kernel (`grri`) + 2. Defining a grid of points (`R_grid`, `Z_grid`) where the field is to be calculated + 3. Calling `_pickup_field` to compute the magnetic field components on that grid using the kernel + and the source perturbation `Bn` # Arguments -- `inputs::VacuumInput`: Struct containing vacuum calculation parameters (n, mpert, mtheta, etc.) -- `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry and basis functions -- `wall::WallGeometry`: Struct with wall geometry -- `Bn::Vector{<:Number}`: Complex vector of Fourier harmonics of the normal magnetic field - perturbation at the plasma surface, `B_n = B_n_real + i*B_n_imag`. Length must be `mpert`. -- `R_grid::AbstractVector`: Vector of R coordinates for the output field grid -- `Z_grid::AbstractVector`: Vector of Z coordinates for the output field grid + - `inputs::VacuumInput`: Struct containing vacuum calculation parameters (n, mpert, mtheta, etc.) + - `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry and basis functions + - `wall::WallGeometry`: Struct with wall geometry + - `Bn::Vector{<:Number}`: Complex vector of Fourier harmonics of the normal magnetic field + perturbation at the plasma surface, `B_n = B_n_real + i*B_n_imag`. Length must be `mpert`. + - `R_grid::AbstractVector`: Vector of R coordinates for the output field grid + - `Z_grid::AbstractVector`: Vector of Z coordinates for the output field grid # Returns -- `B_R::Matrix{ComplexF64}`: The R-component of the magnetic field on the grid -- `B_Z::Matrix{ComplexF64}`: The Z-component of the magnetic field on the grid -- `B_phi::Matrix{ComplexF64}`: The toroidal component of the magnetic field on the grid -- `grid_info::Matrix{Int}`: Information about the grid points (1=inside plasma, 0=outside) + - `B_R::Matrix{ComplexF64}`: The R-component of the magnetic field on the grid + - `B_Z::Matrix{ComplexF64}`: The Z-component of the magnetic field on the grid + - `B_phi::Matrix{ComplexF64}`: The toroidal component of the magnetic field on the grid + - `grid_info::Matrix{Int}`: Information about the grid points (1=inside plasma, 0=outside) """ -function compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall::WallGeometry, - Bn::Vector{<:Number}, R_grid::AbstractVector, Z_grid::AbstractVector) +function compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall::WallGeometry, + Bn::Vector{<:Number}, R_grid::AbstractVector, Z_grid::AbstractVector) # 1. Call vaccal! to get the inverted Green's function matrix - # The Fortran version calls the whole chain (ent33 -> vaccal), + # The Fortran version calls the whole chain (ent33 -> vaccal), # here we assume vaccal! provides what we need. wv, grri = vaccal!(inputs, plasma_surf, wall) @@ -367,7 +370,7 @@ function compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, # 2. Define grid and parameters for pickup routine nx = length(R_grid) nz = length(Z_grid) - + # 3. Call the field pickup routine B_R, B_Z, B_phi, grid_info = _pickup_field( inputs, plasma_surf, grri, Bn_real, Bn_imag, R_grid, Z_grid @@ -375,4 +378,4 @@ function compute_vacuum_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, return B_R, B_Z, B_phi, grid_info end -end \ No newline at end of file +end diff --git a/src/Vacuum/VacuumInternals.jl b/src/Vacuum/VacuumInternals.jl index 3868a7946..12810aa60 100644 --- a/src/Vacuum/VacuumInternals.jl +++ b/src/Vacuum/VacuumInternals.jl @@ -1,9 +1,9 @@ # Gaussian quadrature weights and points for 8-point integration (used for kernel! function) const GAUSSIANWEIGHTS = [0.101228536290376, 0.222381034453374, 0.313706645877887, 0.362683783378362, - 0.362683783378362, 0.313706645877887, 0.222381034453374, 0.101228536290376] + 0.362683783378362, 0.313706645877887, 0.222381034453374, 0.101228536290376] const GAUSSIANPOINTS = [-0.960289856497536, -0.796666477413627, -0.525532409916329, -0.183434642495650, - 0.183434642495650, 0.525532409916329, 0.796666477413627, 0.960289856497536] + 0.183434642495650, 0.525532409916329, 0.796666477413627, 0.960289856497536] # 32-point Gaussian quadrature abscissae (used for Pn_minus_half_2007 function when nρ̂>0.1) const GAUSSIANWEIGHTS32 = [ @@ -34,60 +34,69 @@ const GAUSSIANPOINTS32 = [ -0.506899908932229390024, -0.421351276130635345364, -0.331868602282127649780, -0.239287362252137074545, -0.144471961582796493485, -0.048307665687738316235, - 0.048307665687738316235, 0.144471961582796493485, - 0.239287362252137074545, 0.331868602282127649780, - 0.421351276130635345364, 0.506899908932229390024, - 0.587715757240762329041, 0.663044266930215200975, - 0.732182118740289680387, 0.794483795967942406963, - 0.849367613732569970134, 0.896321155766052123965, - 0.934906075937739689171, 0.964762255587506430774, - 0.985611511545268335400, 0.997263861849481563545 + 0.048307665687738316235, 0.144471961582796493485, + 0.239287362252137074545, 0.331868602282127649780, + 0.421351276130635345364, 0.506899908932229390024, + 0.587715757240762329041, 0.663044266930215200975, + 0.732182118740289680387, 0.794483795967942406963, + 0.849367613732569970134, 0.896321155766052123965, + 0.934906075937739689171, 0.964762255587506430774, + 0.985611511545268335400, 0.997263861849481563545 ] """ - kernel!(grad_greenfunction_mat, greenfunction_mat, x_obspoints, z_obspoints, x_sourcepoints, z_sourcepoints, j1, j2, isgn, iopw, iops, inputs; xwall=nothing, zwall=nothing) + kernel!(grad_greenfunction_mat, greenfunction_mat, x_obspoints, z_obspoints, x_sourcepoints, z_sourcepoints, j1, j2, isgn, iops, inputs; xwall=nothing, zwall=nothing) Compute kernels of integral equation for Laplace's equation in a torus. -**WARNING: This kernel only supports closed toroidal walls currently. +**WARNING: This kernel only supports closed toroidal walls currently. The residue calculation needs to be updated for open walls.** # Arguments -- `grad_greenfunction_mat`: Gradient Green's function matrix (output) -- `greenfunction_mat`: Green's function matrix (output) -- `x_obspoints`: Observer x coordinates (R coordinates) -- `z_obspoints`: Observer z coordinates (Z coordinates) -- `x_sourcepoints`: Source x coordinates (R coordinates) -- `z_sourcepoints`: Source z coordinates (Z coordinates) -- `j1`: Boundary condition index for observer (1=plasma, 2=wall) -- `j2`: Boundary condition index for source (1=plasma, 2=wall) -- `iopw`: Wall option (0=inactive, 1=active) -- `n`: Toroidal mode number - + - `grad_greenfunction_mat`: Gradient Green's function matrix (output) + - `greenfunction_mat`: Green's function matrix (output) + - `x_obspoints`: Observer x coordinates (R coordinates) + - `z_obspoints`: Observer z coordinates (Z coordinates) + - `x_sourcepoints`: Source x coordinates (R coordinates) + - `z_sourcepoints`: Source z coordinates (Z coordinates) + - `j1/j2`: Block index for observer/source (1=plasma, 2=wall) + - `n`: Toroidal mode number # Returns Modifies `grad_greenfunction_mat` and `greenfunction_mat` in place. +Note that greenfunction_mat is zeroed each time this function is called, +but grad_greenfunction_mat is not since it fills a different block of the +(2 * mtheta, 2 * mtheta) depending on the source/observer. # Notes -- Uses Simpson's rule for integration away from singular points -- Uses Gaussian quadrature near singular points for improved accuracy -- Implements analytical singularity removal following Chance 1997 + - Uses Simpson's rule for integration away from singular points + - Uses Gaussian quadrature near singular points for improved accuracy + - Implements analytical singularity removal following Chance 1997 """ -function kernel!(grad_greenfunction_mat::Matrix{Float64}, greenfunction_mat::Matrix{Float64}, x_obspoints::Vector{Float64}, z_obspoints::Vector{Float64}, x_sourcepoints::Vector{Float64}, z_sourcepoints::Vector{Float64}, j1::Int, j2::Int, iopw::Int, n::Int) +function kernel!( + grad_greenfunction_mat::Matrix{Float64}, + greenfunction_mat::Matrix{Float64}, + x_obspoints::Vector{Float64}, + z_obspoints::Vector{Float64}, + x_sourcepoints::Vector{Float64}, + z_sourcepoints::Vector{Float64}, + j1::Int, + j2::Int, + n::Int +) # These used to be function arguments, but can just set inside here based on j1/j2 - isgn = 2 * j2 - 3 - iops = (j1 == 1 && j2 == 1) ? 1 : 0 # only equal to 1 for plasma/plasma block + plasma_plasma_block = j1 == 1 && j2 == 1 # previously iops + plasma_is_source = j2 == 1 # previously iopw + isgn = plasma_is_source ? -1 : 1 mtheta = length(x_obspoints) dtheta = 2π / mtheta - ak0i = 0.0 - jres = 1 - theta_grid = range(start=0; length=mtheta, step=dtheta) + theta_grid = range(; start=0, length=mtheta, step=dtheta) # Zero out greenfunction_mat at start of each kernel call (matches Fortran behavior) fill!(greenfunction_mat, 0.0) @@ -96,72 +105,55 @@ function kernel!(grad_greenfunction_mat::Matrix{Float64}, greenfunction_mat::Mat error("Length of input arrays (xobs, zobs, xsource, zsce) are different. All length should be the same") end - # definition for solving parameters - simpson_b1=1*dtheta/3 - simpson_b2=2*dtheta/3 - simpson_b4=4*dtheta/3 - - simpson_a1=1*dtheta/3 - simpson_a2=2*dtheta/3 - simpson_a4=4*dtheta/3 - - slog1m = 3*dtheta*(log(dtheta)-1/3) - slog1p = slog1m - - log2dtheta = log(2*dtheta) - log_correction_0=16.0*dtheta*(log2dtheta-68.0/15.0)/15.0 - log_correction_1=128.0*dtheta*(log2dtheta-8.0/15.0)/45.0 - log_correction_2=4.0*dtheta*(7.0*log2dtheta-11.0/15.0)/45.0 + # S₁ᵢ in Chance 1997, eq.(78) + log_correction_0=16.0*dtheta*(log(2*dtheta)-68.0/15.0)/15.0 + log_correction_1=128.0*dtheta*(log(2*dtheta)-8.0/15.0)/45.0 + log_correction_2=4.0*dtheta*(7.0*log(2*dtheta)-11.0/15.0)/45.0 has_zero_crossing = false - # See if the conductor surface needs to be checked for singular points - if j1 == 2 || j2 == 2 + # See if the conductor surface crosses R=0 (x=0) anywhere + if !plasma_plasma_block # Initialize jbot and jtop, and figure out wall geometry based on which block we are in jbot = jtop = mtheta/2+1 - xwall = j2 == 2 ? x_sourcepoints : x_obspoints - - # Find where sign of wall x point acrosses zero. - # has_zero_crossing means there is 0-crossing point + xwall = plasma_is_source ? x_obspoints : x_sourcepoints + # Find if sign of wall x point crosses zero for i in 1:mtheta next_i = i == mtheta ? 1 : i + 1 if xwall[i] * xwall[next_i] <= 0.0 - jbot = xwall[i] > 0.0 ? i : jbot - jtop = xwall[i] < 0.0 ? i + 1 : jtop + jbot = xwall[i] > 0.0 ? i : # index where the wall leaves positive R + jtop = xwall[i] < 0.0 ? i + 1 : jtop # index where the wall returns to positive R has_zero_crossing = true end end end - # do spline and calc derivative for Z'_θ and X'_θ in eq.(51) - spline_x = cubic_spline_interpolation(theta_grid, x_sourcepoints, extrapolation_bc=Interpolations.Periodic()) - spline_z = cubic_spline_interpolation(theta_grid, z_sourcepoints, extrapolation_bc=Interpolations.Periodic()) - - gradients_x = (t -> Interpolations.gradient(spline_x, t)).(theta_grid) - gradients_z = (t -> Interpolations.gradient(spline_z, t)).(theta_grid) - dx_dtheta = first.(gradients_x) # d x / d theta - dz_dtheta = first.(gradients_z) # d z / d theta + # Used for Z'_θ and X'_θ in eq.(51) + spline_x = cubic_spline_interpolation(theta_grid, x_sourcepoints; extrapolation_bc=Interpolations.Periodic()) + spline_z = cubic_spline_interpolation(theta_grid, z_sourcepoints; extrapolation_bc=Interpolations.Periodic()) + dx_dtheta = [Interpolations.gradient(spline_x, t)[1] for t in theta_grid] + dz_dtheta = [Interpolations.gradient(spline_z, t)[1] for t in theta_grid] - # observer loop - work = zeros(mtheta) + # Loop through observer points for j in 1:mtheta - # initialize variable - x_obs=x_obspoints[j] #observation point + # Initialize variables + x_obs=x_obspoints[j] z_obs=z_obspoints[j] - theta_obs=theta_grid[j] # theta value - fill!(work, 0.0) - - # if the point of observation point is in negative, We cannot use green func - # This is same for source point + theta_obs=theta_grid[j] + grad_green_0 = 0.0 # simpson integral for coupling_0 (𝒥 ∇'𝒢⁰∇'ℒ) + # Take a view of the appropriate row of grad_greenfunction_mat for this observer point + grad_green_work = @view(grad_greenfunction_mat[(j1-1)*mtheta+j, (j2-1)*mtheta .+ (1:mtheta)]) + + # if observation point is negative, we cannot use green function if x_obs < 0.0 - (j2 == 2) && (work[j] = 1.0) + !plasma_is_source && (grad_green_work[j] = 1.0) continue end - - # set istart and iend - iend = 2 # end point of integration - aval1=0.0 # ∇𝒢_0 - # if wall crossing zero and wall is source, iend set. - if has_zero_crossing && j2 == 2 + + # Compute istart and iend (start/end index of integration to avoid singularity) + # If no zero crossing, istart = iend = 2 + iend = 2 + if has_zero_crossing && !plasma_is_source + # Determine iend/istart so Simpson sweep avoids integrating across the (x=0) discontinuity near the observer index j if jbot - j == 1 iend = 3 elseif jbot - j == 0 @@ -172,167 +164,117 @@ function kernel!(grad_greenfunction_mat::Matrix{Float64}, greenfunction_mat::Mat iend = 1 end end - istart = 4 - iend # starting point of integration - - # 4-3 on mth. then mths is equals to mtheta+3 ??? I'm not sure - mth_source=mtheta-(istart+iend-1) + istart = 4 - iend - # loop for source point - for i in 1:mth_source - - # 4.5 get source point index(ic) theta(theta), X(xt), and Z(zt) - ic = i + j + istart - 1 + # Perform Simpson integration for nonsingular source points (excludes 3 singular points) + # For cases where wall doesn't cross x=0 (iend = istart = 2), the singular points are j-1, j, j+1 + for ic in (j+istart):(j+mtheta-iend) + # Wrap source point index (ic) so it ranges from [1, mtheta] if ic > mtheta ic = ic - mtheta end - # theta_source=(ic-1)*dtheta x_source=x_sourcepoints[ic] z_source=z_sourcepoints[ic] - # if source point is in negative, we cannot use green function + # if source point is negative, we cannot use green function # & if source point(ic) and obs point (j) is same, it's singular (x_source < 0 || ic == j) && continue - # Call green function # G_n is 2pi𝒢ⁿ; coupling_n is 𝒥 ∇'𝒢ⁿ∇'ℒ; coupling_0 is 𝒥 ∇'𝒢ⁿ∇'ℒ for n=0 - G_n, coupling_n, coupling_0 = green(x_obs,z_obs,x_source,z_source, - dx_dtheta[ic],dz_dtheta[ic],n) + G_n, coupling_n, coupling_0 = green(x_obs, z_obs, x_source, z_source, dx_dtheta[ic], dz_dtheta[ic], n) - # simpson integral. 4 for odd, 2 for even, and 1 for others. - simpson_b=simpson_b2 - if (i÷2)*2 == i - simpson_b=simpson_b4 - end - if (i == 1)||(i == mth_source) - simpson_b=simpson_b1 - end - simpson_a=simpson_a2 - if (i÷2)*2 == i - simpson_a=simpson_a4 - end - if (i == 1)||(i == mth_source) - simpson_a=simpson_a1 - end + # Compute composite Simpson's 1/3 rule weight (https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_1/3_rule) + endpoint = (ic == j + istart || ic == j + mtheta - iend) + odd_index = isodd(ic - j - istart) + wsimpson = (endpoint ? 1 : (odd_index ? 4 : 2)) * dtheta / 3 - # work and gren - # work : simpson integral for coupling_n (𝒥 ∇'𝒢ⁿ∇'ℒ) - # gren : log singularity values accumulated. simpson integral for G_n - # aval1 : aval1 - work[ic] += isgn*coupling_n*simpson_a - greenfunction_mat[j,ic] += G_n*simpson_b # integral of G_n (log singularity) - aval1 += coupling_0 * simpson_a + # Sum contributions to Green's function matrices using Simpson weight + grad_green_work[ic] += isgn * coupling_n * wsimpson + greenfunction_mat[j, ic] += G_n * wsimpson + grad_green_0 += coupling_0 * wsimpson end - - # if it's plasma/plasma, wall/wall and in negative wall point, skip loop - # obs : j1 = 1(plasma), wall(2) - # src : j1 = 1(plasma), wall(2) - if (j1+j2) != 2 && has_zero_crossing && j > jbot && j < jtop + + # Skip singularity calculation for R < 0 wall points + if has_zero_crossing && j > jbot && j < jtop continue end - # Get js - js1=mod(j-iend+mtheta-1,mtheta)+1 - js2=mod(j-iend+mtheta,mtheta)+1 - js3=mod(j-iend+mtheta+1,mtheta)+1 - js4=mod(j-iend+mtheta+2,mtheta)+1 - js5=mod(j-iend+mtheta+3,mtheta)+1 - - # Singular points when source point and obs point are the same - # integration each left and right by gaussian quadrature - for ilr in [1,2] + # Perform Gaussian quadrature for singular points (source = obs point) + # Get indices of the singularity region ([j-2, j-1, j, j+1, j+2] for iend = 2) + js = mod.(j - iend .+ ((mtheta-1):(mtheta+3)), mtheta) .+ 1 + # Integrate region of length 2 * dtheta on left (ilr = 1)/right (ilr = 2) of singularity + for ilr in [1, 2] gauss_xleft = theta_obs + (2*ilr-iend-2)*dtheta gauss_xright = gauss_xleft + 2 * dtheta gauss_xavg = (gauss_xright + gauss_xleft)/2 - gauss_halfdifference = (gauss_xright - gauss_xleft)/2 - tgaus = gauss_xavg .+ GAUSSIANPOINTS .* gauss_halfdifference # tgaus is 8 point gauss points, since GAUSSIANPOINTS is for only [-1,1] - # for each of 8 gaussian points - for ig in 1:8 - tgaus0 = tgaus[ig] #i-th value of for 8 points, in theta - tgaus0 = mod(tgaus0, 2π) - - # get X, X', Z, Z' for gaussian point - xt = spline_x(tgaus0) - xtp = Interpolations.gradient(spline_x, tgaus0)[1] - zt = spline_z(tgaus0) - ztp = Interpolations.gradient(spline_z, tgaus0)[1] - - # call green function - G_n, coupling_n, coupling_0 = green(x_obs,z_obs,xt,zt,xtp,ztp,n) - - # add logarithm to G_n to analytically isolate the singularity. Chance eq.(75) - # iops = 1 - G_n_nonsingular = G_n + iops * log((theta_obs-tgaus[ig])^2)/x_obs - - # calc GAUSSIANWEIGHTS. GAUSSIANWEIGHTS contains gaussian quadrature weights for each 8 points - wgbg=GAUSSIANWEIGHTS[ig]*gauss_halfdifference - - # calc pgaus. below Chance eq.(77) - # All are the noted index, times (wgbg)/(2 pgaus) - # A0 has no extra factor - # A1 has an extra factor of (pgaus +/- 1) - # A2 has an extra factor of (pgaus +/- 2) - pgaus=(tgaus[ig]-theta_obs-(2-iend)*dtheta)/dtheta - pgaus2=pgaus*pgaus - A0 = (pgaus2-1)*(pgaus2-4)/4.0 * wgbg - A1_plus = -(pgaus+1)*pgaus*(pgaus2-4)/6.0 * wgbg - A1_minus = -(pgaus-1)*pgaus*(pgaus2-4)/6.0 * wgbg - A2_plus = (pgaus2-1)*pgaus*(pgaus+2)/24.0 * wgbg - A2_minus = (pgaus2-1)*pgaus*(pgaus-2)/24.0 * wgbg - - # add up in work - work[js1] += isgn * coupling_n * A2_minus - work[js2] += isgn * coupling_n * A1_minus - work[js3] += isgn * coupling_n * A0 - work[js4] += isgn * coupling_n * A1_plus - work[js5] += isgn * coupling_n * A2_plus - - # minus diverging value - work[j] -= isgn * coupling_0 * wgbg - if j == jres - ak0i -= isgn * coupling_0 * wgbg + theta_gauss = gauss_xavg .+ GAUSSIANPOINTS .* dtheta # tgaus is 8 point gauss points, since GAUSSIANPOINTS is for only [-1,1] + for ig in 1:8 # 8-point Gaussian quadrature + # Compute green function for this Gaussian point + theta_gauss0 = mod(theta_gauss[ig], 2π) + x_gauss = spline_x(theta_gauss0) + dx_dtheta_gauss = Interpolations.gradient(spline_x, theta_gauss0)[1] + z_gauss = spline_z(theta_gauss0) + dz_dtheta_gauss = Interpolations.gradient(spline_z, theta_gauss0)[1] + G_n, coupling_n, coupling_0 = green(x_obs, z_obs, x_gauss, z_gauss, dx_dtheta_gauss, dz_dtheta_gauss, n) + + # Add logarithm to G_n to analytically isolate the singularity (first type), Chance eq.(75) + G_n_nonsingular = plasma_plasma_block ? G_n + log((theta_obs-theta_gauss[ig])^2)/x_obs : G_n + + # Redefine hardcoded Gaussian weights on the interval [-1, 1] to physical interval with length 2 * dtheta + wgauss = GAUSSIANWEIGHTS[ig] * dtheta + # Calculate p = θ/Δ = (θⱼ - θ')/Δ, 0 at observation point, ±1,±2 at other 5-point stencil nodes + pgauss=(theta_gauss[ig]-theta_obs-(2-iend)*dtheta)/dtheta + # Compute 5-point Lagrange basis polynomials at the Gauss point and multiply by quadrature weight + A0 = (pgauss^2-1)*(pgauss^2-4)/4.0 * wgauss + A1_plus = -(pgauss+1)*pgauss*(pgauss^2-4)/6.0 * wgauss + A1_minus = -(pgauss-1)*pgauss*(pgauss^2-4)/6.0 * wgauss + A2_plus = (pgauss^2-1)*pgauss*(pgauss+2)/24.0 * wgauss + A2_minus = (pgauss^2-1)*pgauss*(pgauss-2)/24.0 * wgauss + + # First type of singularity: 𝒢ⁿ, occurs plasma as source only (see RHS of Chance eqs. 26/27) + if plasma_is_source + greenfunction_mat[j, js[1]] += G_n_nonsingular * A2_minus + greenfunction_mat[j, js[2]] += G_n_nonsingular * A1_minus + greenfunction_mat[j, js[3]] += G_n_nonsingular * A0 + greenfunction_mat[j, js[4]] += G_n_nonsingular * A1_plus + greenfunction_mat[j, js[5]] += G_n_nonsingular * A2_plus end - - # skip when plasma, no skip when considering wall - (iopw == 0) && continue - - # if Wall is considered(wall/wall, plasma/wall, wall/plasma), add up G_n_nonsingular value - greenfunction_mat[j,js1] += G_n_nonsingular * A2_minus - greenfunction_mat[j,js2] += G_n_nonsingular * A1_minus - greenfunction_mat[j,js3] += G_n_nonsingular * A0 - greenfunction_mat[j,js4] += G_n_nonsingular * A1_plus - greenfunction_mat[j,js5] += G_n_nonsingular * A2_plus + + # Second type of singularity: 𝒦ⁿ + # Eq. 86: 𝒦ⁿαᵢ - δⱼᵢK⁰ (js[3] = j if iend=2) + grad_green_work[js[1]] += isgn * coupling_n * A2_minus + grad_green_work[js[2]] += isgn * coupling_n * A1_minus + grad_green_work[js[3]] += isgn * coupling_n * A0 + grad_green_work[js[4]] += isgn * coupling_n * A1_plus + grad_green_work[js[5]] += isgn * coupling_n * A2_plus + # Subtract off the diverging singular n=0 component + grad_green_work[j] -= isgn * coupling_0 * wgauss end end - # Set residue based on logic similar to Table I of Chance 1997 - residue = (j1 == 2) ? 2.0 : 0.0 + # Set residue based on logic similar to Table I of Chance 1997 + existing δⱼᵢ in eq. 69 # Would need to pass in wall geometry to generalize this to open walls is_closed_toroidal = true if is_closed_toroidal - residue_diagonal=(2-j1)*(2-j2)+(j1-1)*(j2-1) - residue_k0=(2-j1)*(2-j2)+(j1-3)*(j2-1) - residue=residue_diagonal+residue_k0 - end - - # minus residue value - work[j] = work[j] - isgn * aval1 + residue - if j == jres - ak0i -= isgn * aval1 + residue = (j1 == 2.0) ? 0.0 : (j2 == 1 ? 2.0 : -2.0) # Chance eq. 89 + else + # TODO: this line can be gotten rid of if we are never doing open walls + residue = (j1 == j2) ? 2.0 : 0.0 # Chance eq. 90 end - - # Only when plasma/plasma, log singularity activate. (S1) - if iops == 1 && iopw != 0 - greenfunction_mat[j,js1] -= log_correction_2 / x_obs - greenfunction_mat[j,js2] -= log_correction_1 / x_obs - greenfunction_mat[j,js3] -= log_correction_0 / x_obs - greenfunction_mat[j,js4] -= log_correction_1/x_obs - greenfunction_mat[j,js5] -= log_correction_2/x_obs + # Subtract regular integral component of δⱼᵢK⁰ in eq. 83 and add residue value in eq. 89/90 + grad_green_work[j] = grad_green_work[j] - isgn * grad_green_0 + residue + + # Subtract off analytic singular integral from Chance eq.(75) if plasma-plasma block + if plasma_plasma_block + greenfunction_mat[j, js[1]] -= log_correction_2 / x_obs + greenfunction_mat[j, js[2]] -= log_correction_1 / x_obs + greenfunction_mat[j, js[3]] -= log_correction_0 / x_obs + greenfunction_mat[j, js[4]] -= log_correction_1 / x_obs + greenfunction_mat[j, js[5]] -= log_correction_2 / x_obs end - - # Store all the datas of work in grdgre, gren - @views grad_greenfunction_mat[(j1-1)*mtheta + j, (j2-1)*mtheta .+ (1:mtheta)] .= work - @views greenfunction_mat[j, 1:mtheta] ./= 2π end + # Since we computed 2π𝒢, divide by 2π to get 𝒢 + greenfunction_mat ./= 2π end """ @@ -342,20 +284,20 @@ Perform the inverse Fourier transform of `gil` onto `gll` using Fourier coeffici # Arguments -- `gll`: Output matrix (mpert × mpert) updated in-place -- `gil`: Input matrix (mtheta × mpert) containing Fourier-space data -- `cs`: Fourier coefficient matrix (mtheta × mpert) -- `m00`: Integer offset in the gil matrix (row offset) -- `l00`: Integer offset in the gil matrix (column offset) + - `gll`: Output matrix (mpert × mpert) updated in-place + - `gil`: Input matrix (mtheta × mpert) containing Fourier-space data + - `cs`: Fourier coefficient matrix (mtheta × mpert) + - `m00`: Integer offset in the gil matrix (row offset) + - `l00`: Integer offset in the gil matrix (column offset) # Notes -- Computes: `gll[l2, l1] = (2π * dth) * Σ_i cs[i, l2] * gil[i, l1]` -- Performs the same function as fouranv in the Fortran code. + - Computes: `gll[l2, l1] = (2π * dth) * Σ_i cs[i, l2] * gil[i, l1]` + - Performs the same function as fouranv in the Fortran code. # Returns -- gll(l2,l1) : output matrix updated in-place (mpert × mpert) + - gll(l2,l1) : output matrix updated in-place (mpert × mpert) """ function fourier_inverse_transform!(gll::Matrix{Float64}, gil::Matrix{Float64}, cs::Matrix{Float64}, m00::Int, l00::Int) @@ -366,7 +308,7 @@ function fourier_inverse_transform!(gll::Matrix{Float64}, gil::Matrix{Float64}, # Inverse Fourier transform via matrix multiply: gll = cs^T * gil * (2π * dth) # This computes: gll[l2, l1] = (2π * dth) * Σ_i cs[i, l2] * gil[i, l1] dth = 2π / mtheta - mul!(gll, cs', view(gil, m00+1:m00+mtheta, l00+1:l00+mpert), 2π * dth, 0.0) + mul!(gll, cs', view(gil, (m00+1):(m00+mtheta), (l00+1):(l00+mpert)), 2π * dth, 0.0) end """ @@ -375,7 +317,7 @@ end Purpose: This routine performs a truncated Fourier transform of gij onto gil using Fourier coefficients stored in cs. - + Inputs: gij(i,j) : input matrix of size (mth × mth), the "physical-space" data cs(j,l) : Fourier coefficient matrix (mth × mpert) @@ -390,10 +332,10 @@ function fourier_transform!(gil::Matrix{Float64}, gij::Matrix{Float64}, cs::Matr # Zero out relevant gil block mtheta, mpert = size(cs) - fill!(view(gil, m00+1:m00+mtheta, l00+1:l00+mpert), 0.0) + fill!(view(gil, (m00+1):(m00+mtheta), (l00+1):(l00+mpert)), 0.0) # Fourier transform via matrix multiply: gil[i, l] = Σ_j gij[i, j] * cs[j, l] - mul!(view(gil, m00+1:m00+mtheta, l00+1:l00+mpert), gij, cs) + mul!(view(gil, (m00+1):(m00+mtheta), (l00+1):(l00+mpert)), gij, cs) end # Returns the array of derivatives at all x points, I think this acts like difspl @@ -411,22 +353,22 @@ Replaces `lagp`, `lagpe4`, and `lag` from Fortran code. # Arguments -- `ax::Vector{Float64}`: Array of x-coordinates for the interpolation points -- `af::Vector{Float64}`: Array of y-coordinates (function values) for the interpolation points -- `m::Int`: Number of interpolation points -- `nl::Int`: Number of points to use for the local interpolation (polynomial degree + 1) -- `x::Float64`: The x-value at which to evaluate the interpolated function and/or its derivative -- `iop::Int`: Flag controlling output (0 = value only, 1 = value and derivative) + - `ax::Vector{Float64}`: Array of x-coordinates for the interpolation points + - `af::Vector{Float64}`: Array of y-coordinates (function values) for the interpolation points + - `m::Int`: Number of interpolation points + - `nl::Int`: Number of points to use for the local interpolation (polynomial degree + 1) + - `x::Float64`: The x-value at which to evaluate the interpolated function and/or its derivative + - `iop::Int`: Flag controlling output (0 = value only, 1 = value and derivative) # Returns -- `f::Float64`: The interpolated function value at `x` -- `df::Float64`: The interpolated function derivative at `x` (0.0 if iop=0) + - `f::Float64`: The interpolated function value at `x` + - `df::Float64`: The interpolated function derivative at `x` (0.0 if iop=0) # Notes -- Uses local Lagrange interpolation with `nl` points centered around `x` -- Automatically adjusts interpolation window to stay within array bounds + - Uses local Lagrange interpolation with `nl` points centered around `x` + - Automatically adjusts interpolation window to stay within array bounds """ function lagrange1d(ax::Vector{Float64}, af::Vector{Float64}, m::Int, nl::Int, x::Float64, iop::Int) @@ -437,7 +379,7 @@ function lagrange1d(ax::Vector{Float64}, af::Vector{Float64}, m::Int, nl::Int, x jn = findfirst(i -> ax[i] >= x, 1:m) jn = (jn === nothing) ? m : jn jn = max(jn - 1, 1) - if jn < m && abs(ax[jn + 1] - x) < abs(x - ax[jn]) + if jn < m && abs(ax[jn+1] - x) < abs(x - ax[jn]) jn += 1 end @@ -499,25 +441,25 @@ end Resample the input array `vecin` using a periodic cubic spline to an output array of length `mtheta`. -This function unifies the Fortran functions `trans`, `transdx`, and `transdxx` into a single +This function unifies the Fortran functions `trans`, `transdx`, and `transdxx` into a single function with optional offset parameters. # Arguments -- `vecin::Vector{Float64}`: Input array to be resampled -- `mtheta::Int`: Desired length of the output array -- `dx0::Float64`: Global offset added to all x-coordinates (default 0, applied as `x += dx0 / mtheta_in`) -- `dx1::Float64`: Fine offset added to each index (default 0, applied as `ai = (i-1) + dx1`) + - `vecin::Vector{Float64}`: Input array to be resampled + - `mtheta::Int`: Desired length of the output array + - `dx0::Float64`: Global offset added to all x-coordinates (default 0, applied as `x += dx0 / mtheta_in`) + - `dx1::Float64`: Fine offset added to each index (default 0, applied as `ai = (i-1) + dx1`) # Returns -- `vecout::Vector{Float64}`: The resampled output array (length `mtheta`) + - `vecout::Vector{Float64}`: The resampled output array (length `mtheta`) # Notes -- If `mtheta == length(vecin)`, returns the input vector unchanged -- Uses periodic cubic spline interpolation for resampling -- Input grid is normalized to [0, 1] for interpolation + - If `mtheta == length(vecin)`, returns the input vector unchanged + - Uses periodic cubic spline interpolation for resampling + - Input grid is normalized to [0, 1] for interpolation """ function interp_to_new_grid(vecin::Vector{Float64}, mtheta::Int; dx0=0.0, dx1=0.0) @@ -627,25 +569,25 @@ end """ Pn_minus_half_1997(s, n) -Compute the Legendre function of the first kind of order -1/2, P^n_{-1/2}(s), +Compute the Legendre function of the first kind of order -1/2, P^n_{-1/2}(s), recursively using Chance 1997 equations (47)-(50). -The implementation follows the original Fortran code. Note: equation (50) in the paper +The implementation follows the original Fortran code. Note: equation (50) in the paper has a typo where the exponent should be -1/4 instead of +1/2. # Arguments -- `s::Real`: Legendre function parameter (s > 1) -- `n::Int`: Maximum order n (n ≥ 0) + - `s::Real`: Legendre function parameter (s > 1) + - `n::Int`: Maximum order n (n ≥ 0) # Returns -- `P::Vector{Float64}`: Array of values P^0_{-1/2}(s) through P^{n+1}_{-1/2}(s) + - `P::Vector{Float64}`: Array of values P^0_{-1/2}(s) through P^{n+1}_{-1/2}(s) # Notes -- Uses recursive relation from Chance 1997 eq. (47) -- Base cases computed from eqs. (48)-(50) using elliptic integrals + - Uses recursive relation from Chance 1997 eq. (47) + - Base cases computed from eqs. (48)-(50) using elliptic integrals """ function Pn_minus_half_1997(s::Real, n::Int) @@ -675,20 +617,23 @@ Compute complete elliptic integrals K(m1) and E(m1) using Bulirsch's algorithm. This is the Julia equivalent of the Fortran `ek3` subroutine. # Arguments -- `m1::Float64`: Complementary parameter (1 - k²), where k is the elliptic modulus -- `error::Float64`: Convergence tolerance (default 1e-8) -- `maxit::Int`: Maximum iterations (default 10) + + - `m1::Float64`: Complementary parameter (1 - k²), where k is the elliptic modulus + - `error::Float64`: Convergence tolerance (default 1e-8) + - `maxit::Int`: Maximum iterations (default 10) # Returns -- `K::Float64`: Complete elliptic integral of the first kind K(m1) -- `E::Float64`: Complete elliptic integral of the second kind E(m1) -- `convergence::Float64`: Convergence metric -- `iterations::Int`: Number of iterations performed + + - `K::Float64`: Complete elliptic integral of the first kind K(m1) + - `E::Float64`: Complete elliptic integral of the second kind E(m1) + - `convergence::Float64`: Convergence metric + - `iterations::Int`: Number of iterations performed # Notes -- Based on Bulirsch's method as described in Numerical Recipes -- Precision is approximately error² -- Reference: JCP 221 (2007) 330-348 + + - Based on Bulirsch's method as described in Numerical Recipes + - Precision is approximately error² + - Reference: JCP 221 (2007) 330-348 """ function elliptic_integrals_bulirsch(m1::Float64; error::Float64=1e-8, maxit::Int=10) @@ -783,21 +728,25 @@ Compute the Legendre function of the first kind of order -1/2, P^n_{-1/2}(s), using methods from Chance J. Comp. Phys 221 (2007) 330-348. This implementation uses: -1. Bulirsch's algorithm for elliptic integrals (more accurate than polynomial approximations) -2. Gaussian integration for large mode numbers (n*rhohat >= 0.1) where rhohat = 1/√(2*y*w) -3. Upward recurrence for small mode numbers + + 1. Bulirsch's algorithm for elliptic integrals (more accurate than polynomial approximations) + 2. Gaussian integration for large mode numbers (n*rhohat >= 0.1) where rhohat = 1/√(2*y*w) + 3. Upward recurrence for small mode numbers # Arguments -- `s::Real`: Legendre function parameter (s > 1) -- `n::Int`: Maximum order n (n ≥ 0) + + - `s::Real`: Legendre function parameter (s > 1) + - `n::Int`: Maximum order n (n ≥ 0) # Returns -- `P::Vector{Float64}`: Array of values P^0_{-1/2}(s) through P^{n+1}_{-1/2}(s) + + - `P::Vector{Float64}`: Array of values P^0_{-1/2}(s) through P^{n+1}_{-1/2}(s) # Notes -- This version is more accurate than Pn_minus_half_1997 for large n -- Expected to diverge from 1997 version at large nloc -- Reference: JCP 221 (2007) 330-348 + + - This version is more accurate than Pn_minus_half_1997 for large n + - Expected to diverge from 1997 version at large nloc + - Reference: JCP 221 (2007) 330-348 # Constants """ function Pn_minus_half_2007(s::Real, n::Int) @@ -826,7 +775,7 @@ function Pn_minus_half_2007(s::Real, n::Int) m1sqrti = sqrt(w) # m1^(-1/4) # Compute elliptic integrals using Bulirsch algorithm - K, E, conv, iters = elliptic_integrals_bulirsch(m1sq, error=1e-15, maxit=20) + K, E, conv, iters = elliptic_integrals_bulirsch(m1sq; error=1e-15, maxit=20) # Base cases: P^0 and P^1 pn = pii * m1sqrt * K @@ -920,27 +869,27 @@ according to equations (36)-(42) of Chance 1997. Replaces `green` from Fortran c # Arguments -- `x_obs`: Observation point R-coordinate (Float64) -- `z_obs`: Observation point Z-coordinate (Float64) -- `x_source`: Source point R-coordinate (Float64) -- `z_source`: Source point Z-coordinate (Float64) -- `dx_dtheta`: Derivative ∂R'/∂θ at source point (Float64) -- `dz_dtheta`: Derivative ∂Z'/∂θ at source point (Float64) -- `n`: Toroidal mode number (Int) -- `uselegacygreenfunction::Bool`: Flag to use the 1997 version of the Legendre function (default false, uses 2007 version) + - `x_obs`: Observation point R-coordinate (Float64) + - `z_obs`: Observation point Z-coordinate (Float64) + - `x_source`: Source point R-coordinate (Float64) + - `z_source`: Source point Z-coordinate (Float64) + - `dx_dtheta`: Derivative ∂R'/∂θ at source point (Float64) + - `dz_dtheta`: Derivative ∂Z'/∂θ at source point (Float64) + - `n`: Toroidal mode number (Int) + - `uselegacygreenfunction::Bool`: Flag to use the 1997 version of the Legendre function (default false, uses 2007 version) # Returns -- `G_n`: 2π𝒢ⁿ(θ,θ′) — Green's function value -- `coupling_n`: 𝒥 ∇'𝒢ⁿ∇'ℒ — Coupling term for mode n -- `coupling_0`: 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ — Coupling term for mode 0 + - `G_n`: 2π𝒢ⁿ(θ,θ′) — Green's function value + - `coupling_n`: 𝒥 ∇'𝒢ⁿ∇'ℒ — Coupling term for mode n + - `coupling_0`: 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ — Coupling term for mode 0 # Notes -- Uses Legendre functions P^n_{-1/2}(s) computed via elliptic integrals -- Implements analytical derivatives from Chance 1997 equations -- The coupling terms include the Jacobian factor from the coordinate transformation -- By default uses the 2007 Legendre function implementation (Bulirsch + Gaussian integration) + - Uses Legendre functions P^n_{-1/2}(s) computed via elliptic integrals + - Implements analytical derivatives from Chance 1997 equations + - The coupling terms include the Jacobian factor from the coordinate transformation + - By default uses the 2007 Legendre function implementation (Bulirsch + Gaussian integration) """ function green(x_obs::Float64, z_obs::Float64, x_source::Float64, z_source::Float64, dx_dtheta::Float64, dz_dtheta::Float64, n::Int; uselegacygreenfunction::Bool=false) @@ -955,15 +904,15 @@ function green(x_obs::Float64, z_obs::Float64, x_source::Float64, z_source::Floa # Chance 1997 eq.(41) ℛ = R R4 = ρ2 * (ρ2 + 4 * x_multiple) - R2 = sqrt(R4) + R2 = sqrt(R4) R = sqrt(R2) R5 = R4 * R # Chance 1997 eq.(42) 𝘴 = s s = (x_obs2 + x_source2 + ζ2) / R2 - - # Legendre functions for - # P⁰ = p0, P¹ = p1, Pⁿ = pn, Pⁿ⁺¹ = pnp1 + + # Legendre functions for + # P⁰ = p0, P¹ = p1, Pⁿ = pn, Pⁿ⁺¹ = pnp1 if uselegacygreenfunction legendre = Pn_minus_half_1997(s, n) else @@ -980,26 +929,26 @@ function green(x_obs::Float64, z_obs::Float64, x_source::Float64, z_source::Floa G_n = gg * pn # Chance 1997 eq.(44) (Note this equation in the paper has an erroneous extra factor of 2π) - grad_gg = gg / R4 /2π + grad_gg = gg / R4 / 2π # ∂Gⁿ/∂X' = dG_dX - xterm1 = (n * (x_obs2 + x_source2 + ζ2)*(x_obs2 - x_source2 + ζ2) - x_source2*(x_source2-x_obs2+ζ2)) * pn - xterm2 = (2.0 * x_source * x_obs * (x_obs2-x_source2+ζ2)) * pnp1 + xterm1 = (n * (x_obs2 + x_source2 + ζ2) * (x_obs2 - x_source2 + ζ2) - x_source2*(x_source2-x_obs2+ζ2)) * pn + xterm2 = (2.0 * x_source * x_obs * (x_obs2-x_source2+ζ2)) * pnp1 dG_dX = grad_gg * (xterm1 + xterm2) / x_source - + # ∂Gⁿ/∂Z' = dG_dZ zterm1 = (2.0 * n + 1.0) * (x_obs2 + x_source2 + ζ2) * pn zterm2 = 4.0 * x_multiple * pnp1 dG_dZ = grad_gg * (zterm1 + zterm2) * ζ - - # Chance 1997 eq.(51) + + # Chance 1997 eq.(51) # 𝒥 ∇'𝒢ⁿ∇'ℒ = aval # ∂X'/∂θ = xtp, ∂Z'/∂θ = ztp coupling_n = -x_source * (dz_dtheta * dG_dX - dx_dtheta * dG_dZ) - # for 𝓃⩵0, aval0 = 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ - dG_dX0_R5 = ((2.0 * x_obs * (x_obs2-x_source2+ζ2)) * p1 - x_source*(x_source2-x_obs2+ζ2) * p0 ) - dG_dZ0_R5 = ζ * ((x_obs2 + x_source2 + ζ2) * p0 + 4.0 * x_multiple * p1) + # for 𝓃⩵0, aval0 = 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ + dG_dX0_R5 = ((2.0 * x_obs * (x_obs2-x_source2+ζ2)) * p1 - x_source * (x_source2-x_obs2+ζ2) * p0) + dG_dZ0_R5 = ζ * ((x_obs2 + x_source2 + ζ2) * p0 + 4.0 * x_multiple * p1) coupling_0 = -x_source * (dz_dtheta * dG_dX0_R5 - dx_dtheta * dG_dZ0_R5) / R5 return G_n, coupling_n, coupling_0 end @@ -1012,34 +961,34 @@ end Calculate the magnetic field on a specified grid using finite differencing of the magnetic scalar potential `chi`. -This is the Julia version of the Fortran `pickup` routine. It computes the vacuum +This is the Julia version of the Fortran `pickup` routine. It computes the vacuum magnetic field perturbation at a set of grid points given the plasma surface perturbation. # Arguments -- `inputs::VacuumInput`: Struct containing vacuum calculation parameters -- `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry and basis functions -- `grri::Matrix{Float64}`: Inverted Green's function response matrix from vaccal! -- `Bn_real::Vector{Float64}`: Real part of normal field Fourier harmonics at plasma surface -- `Bn_imag::Vector{Float64}`: Imaginary part of normal field Fourier harmonics at plasma surface -- `R_grid::AbstractVector`: R-coordinates for output field evaluation -- `Z_grid::AbstractVector`: Z-coordinates for output field evaluation + - `inputs::VacuumInput`: Struct containing vacuum calculation parameters + - `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry and basis functions + - `grri::Matrix{Float64}`: Inverted Green's function response matrix from vaccal! + - `Bn_real::Vector{Float64}`: Real part of normal field Fourier harmonics at plasma surface + - `Bn_imag::Vector{Float64}`: Imaginary part of normal field Fourier harmonics at plasma surface + - `R_grid::AbstractVector`: R-coordinates for output field evaluation + - `Z_grid::AbstractVector`: Z-coordinates for output field evaluation # Returns -- `B_R::Matrix{ComplexF64}`: R-component of magnetic field on grid (nx × nz) -- `B_Z::Matrix{ComplexF64}`: Z-component of magnetic field on grid (nx × nz) -- `B_phi::Matrix{ComplexF64}`: Toroidal component of magnetic field on grid (nx × nz) -- `grid_info::Matrix{Int}`: Grid point classification (1=inside plasma, 0=outside) + - `B_R::Matrix{ComplexF64}`: R-component of magnetic field on grid (nx × nz) + - `B_Z::Matrix{ComplexF64}`: Z-component of magnetic field on grid (nx × nz) + - `B_phi::Matrix{ComplexF64}`: Toroidal component of magnetic field on grid (nx × nz) + - `grid_info::Matrix{Int}`: Grid point classification (1=inside plasma, 0=outside) # Notes -- Uses 5-point finite difference stencil for computing field from potential -- Field components computed as: B_R = -∂χ/∂R, B_Z = -∂χ/∂Z, B_φ = inχ/R + - Uses 5-point finite difference stencil for computing field from potential + - Field components computed as: B_R = -∂χ/∂R, B_Z = -∂χ/∂Z, B_φ = inχ/R """ function _pickup_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, grri::Matrix{Float64}, - Bn_real::Vector{Float64}, Bn_imag::Vector{Float64}, - R_grid::AbstractVector, Z_grid::AbstractVector) + Bn_real::Vector{Float64}, Bn_imag::Vector{Float64}, + R_grid::AbstractVector, Z_grid::AbstractVector) nx = length(R_grid) nz = length(Z_grid) @@ -1066,7 +1015,7 @@ function _pickup_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, grri::M Threads.@threads for i in 1:n_points R, Z = R_points[i], Z_points[i] - + # Points for finite difference stencil observe_points = [ (R, Z + del_Z), @@ -1091,15 +1040,15 @@ function _pickup_field(inputs::VacuumInput, plasma_surf::PlasmaGeometry, grri::M # B_R = -d(chi)/dR, B_Z = -d(chi)/dZ br_r = -(chi_r[3, idx] - chi_r[4, idx]) / (2.0 * del_R) br_i = -(chi_i[3, idx] - chi_i[4, idx]) / (2.0 * del_R) - + bz_r = -(chi_r[1, idx] - chi_r[2, idx]) / (2.0 * del_Z) bz_i = -(chi_i[1, idx] - chi_i[2, idx]) / (2.0 * del_Z) - + # B_phi = i*n*chi / R # Bphi = i*n*(chi_r + i*chi_i)/R = (-n*chi_i + i*n*chi_r)/R bphi_r = -inputs.n * chi_i[5, idx] / R_points[idx] - bphi_i = inputs.n * chi_r[5, idx] / R_points[idx] - + bphi_i = inputs.n * chi_r[5, idx] / R_points[idx] + B_R_complex[i, j] = br_r + ifac * br_i B_Z_complex[i, j] = bz_r + ifac * bz_i B_phi_complex[i, j] = bphi_r + ifac * bphi_i @@ -1130,17 +1079,17 @@ This is the Julia version of the Fortran `loops` subroutine. # Arguments -- `R_grid::AbstractVector`: Vector of R-coordinates defining the grid -- `Z_grid::AbstractVector`: Vector of Z-coordinates defining the grid + - `R_grid::AbstractVector`: Vector of R-coordinates defining the grid + - `Z_grid::AbstractVector`: Vector of Z-coordinates defining the grid # Returns -- `R_points::Vector{Float64}`: Flattened array of R-coordinates (length nx*nz) -- `Z_points::Vector{Float64}`: Flattened array of Z-coordinates (length nx*nz) + - `R_points::Vector{Float64}`: Flattened array of R-coordinates (length nx*nz) + - `Z_points::Vector{Float64}`: Flattened array of Z-coordinates (length nx*nz) # Notes -- Grid points are ordered as: [(R[1],Z[1]), (R[1],Z[2]), ..., (R[1],Z[nz]), (R[2],Z[1]), ...] + - Grid points are ordered as: [(R[1],Z[1]), (R[1],Z[2]), ..., (R[1],Z[nz]), (R[2],Z[1]), ...] """ function _create_pickup_grid(R_grid::AbstractVector, Z_grid::AbstractVector) nx = length(R_grid) @@ -1168,33 +1117,33 @@ by integrating the Green's function response with the source perturbation at the # Arguments -- `R_obs::Float64`: R-coordinate of observation point -- `Z_obs::Float64`: Z-coordinate of observation point -- `inputs::VacuumInput`: Struct containing vacuum calculation parameters -- `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry -- `grri::Matrix{Float64}`: Inverted Green's function response matrix -- `Bn_real::Vector{Float64}`: Real part of normal field Fourier harmonics -- `Bn_imag::Vector{Float64}`: Imaginary part of normal field Fourier harmonics + - `R_obs::Float64`: R-coordinate of observation point + - `Z_obs::Float64`: Z-coordinate of observation point + - `inputs::VacuumInput`: Struct containing vacuum calculation parameters + - `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry + - `grri::Matrix{Float64}`: Inverted Green's function response matrix + - `Bn_real::Vector{Float64}`: Real part of normal field Fourier harmonics + - `Bn_imag::Vector{Float64}`: Imaginary part of normal field Fourier harmonics # Returns -- `chi_real::Float64`: Real part of the magnetic scalar potential at (R_obs, Z_obs) -- `chi_imag::Float64`: Imaginary part of the magnetic scalar potential at (R_obs, Z_obs) + - `chi_real::Float64`: Real part of the magnetic scalar potential at (R_obs, Z_obs) + - `chi_imag::Float64`: Imaginary part of the magnetic scalar potential at (R_obs, Z_obs) # Notes -- The potential is computed via Fourier series over poloidal modes -- Includes coupling term from Green's function derivative -- Factor of -0.5 * dtheta applied from Fortran convention + - The potential is computed via Fourier series over poloidal modes + - Includes coupling term from Green's function derivative + - Factor of -0.5 * dtheta applied from Fortran convention """ function _calculate_potential_chi(R_obs::Float64, Z_obs::Float64, - inputs::VacuumInput, plasma_surf::PlasmaGeometry, - grri::Matrix{Float64}, - Bn_real::Vector{Float64}, Bn_imag::Vector{Float64}) - + inputs::VacuumInput, plasma_surf::PlasmaGeometry, + grri::Matrix{Float64}, + Bn_real::Vector{Float64}, Bn_imag::Vector{Float64}) + chi_real = 0.0 chi_imag = 0.0 - + mtheta = inputs.mtheta mpert = inputs.mpert n = inputs.n @@ -1210,8 +1159,8 @@ function _calculate_potential_chi(R_obs::Float64, Z_obs::Float64, for i_theta in 1:mtheta R_src = plasma_surf.x[i_theta] Z_src = plasma_surf.z[i_theta] - - # Call the low-level Green's function calculator. + + # Call the low-level Green's function calculator. # The `green` function returns the Green's function value itself (G_n) and # the coupling terms for mode n and mode 0. G_n, coupling_n, coupling_0 = green(R_obs, Z_obs, R_src, Z_src, plasma_surf.dx_dtheta[i_theta], plasma_surf.dz_dtheta[i_theta], n) @@ -1241,18 +1190,18 @@ function _calculate_potential_chi(R_obs::Float64, Z_obs::Float64, # grri has structure [ (grri_cc, grri_cs), (grri_sc, grri_ss) ] # The indices for chiwc, chiws in Fortran map to columns of grri chi_wc = grri[i_theta, l1] # Real part kernel - chi_ws = grri[i_theta, mpert + l1] # Imaginary part kernel - + chi_ws = grri[i_theta, mpert+l1] # Imaginary part kernel + term_r += g_real[i_theta, l1] * chi_wc - g_imag[i_theta, l1] * chi_ws term_i += g_imag[i_theta, l1] * chi_wc + g_real[i_theta, l1] * chi_ws end - + chi_real += term_r * Bn_real[l1] - term_i * Bn_imag[l1] chi_imag += term_i * Bn_real[l1] + term_r * Bn_imag[l1] end - + # The factor of 0.5 * isg * dth in Fortran # isg is -1 for plasma surface factor = -0.5 * dtheta return chi_real * factor, chi_imag * factor -end \ No newline at end of file +end From d161887cfe04fbc6ccb571b98ff3522a544d727a Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Tue, 6 Jan 2026 13:50:35 -0500 Subject: [PATCH 2/2] VACUUM - BUGFIX - incorrect indexing --- src/Vacuum/Vacuum.jl | 2 -- src/Vacuum/VacuumInternals.jl | 14 +++++++------- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index e4d9da32e..7b77f9bf2 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -313,8 +313,6 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe end wv = complex.(vacmat, vacmti) - display(wv) - # Create xzpts array xzpts = zeros(Float64, inputs.mtheta, 4) @views xzpts[:, 1] .= plasma_surf.x diff --git a/src/Vacuum/VacuumInternals.jl b/src/Vacuum/VacuumInternals.jl index 12810aa60..c1295e1e7 100644 --- a/src/Vacuum/VacuumInternals.jl +++ b/src/Vacuum/VacuumInternals.jl @@ -44,7 +44,6 @@ const GAUSSIANPOINTS32 = [ 0.985611511545268335400, 0.997263861849481563545 ] - """ kernel!(grad_greenfunction_mat, greenfunction_mat, x_obspoints, z_obspoints, x_sourcepoints, z_sourcepoints, j1, j2, isgn, iops, inputs; xwall=nothing, zwall=nothing) @@ -140,7 +139,7 @@ function kernel!( z_obs=z_obspoints[j] theta_obs=theta_grid[j] grad_green_0 = 0.0 # simpson integral for coupling_0 (𝒥 ∇'𝒢⁰∇'ℒ) - # Take a view of the appropriate row of grad_greenfunction_mat for this observer point + # Workspace = view of appropriate row of grad_greenfunction_mat for this observer point grad_green_work = @view(grad_greenfunction_mat[(j1-1)*mtheta+j, (j2-1)*mtheta .+ (1:mtheta)]) # if observation point is negative, we cannot use green function @@ -168,8 +167,9 @@ function kernel!( # Perform Simpson integration for nonsingular source points (excludes 3 singular points) # For cases where wall doesn't cross x=0 (iend = istart = 2), the singular points are j-1, j, j+1 - for ic in (j+istart):(j+mtheta-iend) - # Wrap source point index (ic) so it ranges from [1, mtheta] + for i in 1:(mtheta-3) + # Get source point index (ic) and ensure it is in range [1, mtheta] + ic = i + j + istart - 1 if ic > mtheta ic = ic - mtheta end @@ -184,9 +184,9 @@ function kernel!( G_n, coupling_n, coupling_0 = green(x_obs, z_obs, x_source, z_source, dx_dtheta[ic], dz_dtheta[ic], n) # Compute composite Simpson's 1/3 rule weight (https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_1/3_rule) - endpoint = (ic == j + istart || ic == j + mtheta - iend) - odd_index = isodd(ic - j - istart) - wsimpson = (endpoint ? 1 : (odd_index ? 4 : 2)) * dtheta / 3 + # Note we set to 4 for even/2 for odd since we index from 1 while the formula assumes indexing from 0 + endpoint = (i == 1)||(i == mtheta - 3) + wsimpson = (endpoint ? 1 : (iseven(i) ? 4 : 2)) * dtheta / 3 # Sum contributions to Green's function matrices using Simpson weight grad_green_work[ic] += isgn * coupling_n * wsimpson