Skip to content

GPEC - ENHANCEMENT - div-free option for psi,theta grid fields and ar… - #365

Draft
ysm1887 wants to merge 1 commit into
developfrom
regdiv_arzphifun
Draft

ysm1887 wants to merge 1 commit into
developfrom
regdiv_arzphifun

Conversation

@ysm1887

@ysm1887 ysm1887 commented Nov 2, 2026

Copy link
Copy Markdown
Collaborator

GPEC - ENHANCEMENT - div-free option for psi,theta grid fields and arzphifun grad(zeta) update

For fast-ion calculations a divergence-free perturbed field can matter, so this adds an output option for that and makes the flux-grid vector potential consistent with it. Both changes are output-only (gpout_xbrzphifun, gpout_arzphifun); reg_div_flag defaults to f, with it off the outputs are identical to develop.

1. gpout_arzphifun: grad(zeta) components w(3,1:2)

gpeq.f (l.757-761) computes the poloidal components of grad(zeta) as

w(3,1) = (twopi*r*rfac/jac)*(rzphi%fx(2)*rzphi%fy(3)-rzphi%fx(3)*(1+rzphi%fy(2)))
w(3,2) = (r/(2*rfac*jac))*(rzphi%fx(3)*rzphi%fy(1)-rzphi%fx(1)*rzphi%fy(3))

gpout_arzphifun used w(3,1)=2*rfac*dphi_psi, w(3,2)=dphi_theta/(2*pi*rfac) instead. This PR implements the same quantity there through the w(1,:), w(2,:) already in the routine (sign opposite to gpeq.f because arzphifun applies it with a minus):

w(3,1) = (rzphi%fx(3)*w(1,1) + rzphi%fy(3)*w(2,1))/twopi
w(3,2) = (rzphi%fx(3)*w(1,2) + rzphi%fy(3)*w(2,2))/twopi

Check (KSTAR n=1, psi_N < 0.4): curl(a) vs xbrzphifun b_R / b_z = 168 % / 46 % before, 0.25 % / 0.14 % after; a_phi unchanged.

2. reg_div_flag: the regularized field is not divergence-free

Issue: with reg_flag=t the output field (b_R, b_z, b_phi) does not satisfy div b = 0. The regularization damps bmt and bmz but not bwp, and div b = 0 ties m*bwt - n*bwz to d(bwp)/dpsi. On the psi,theta grid the b_phi implied by div b = 0 from (b_R, b_z) differs from the output b_phi by ~15 % pointwise near the rational surfaces (4 % numerical floor with reg_flag=f).

reg_div_flag=t (default f) lets the user choose, at output time in gpout_xbrzphifun only, a regularization that keeps div b = 0: only the (n,m) projection of (bwt, bwz) is damped,

T   = (n*bwt + m*bwz)/(n^2 + m^2)
bmt = bwt - (1-f)*n*T,   bmz = bwz - (1-f)*m*T

so m*bmt - n*bmz is unchanged. bwp, xi and all other outputs are untouched. Check: pointwise mismatch 0.147 -> 0.044 (= reg-off floor), xi identical, fields finite in the bands.

Files

gpec/gpout.f (fix + flag block), gpec/gpglobal.F, gpec/gpec.f (flag), input/gpec.in. Built from install/ with ifort; added fixed-format lines <= 72 columns. The same approximation also appears in gpdiag_singcurs (gpdiag.f l.476, 502) — not touched here.

…zphifun grad(zeta) update

For fast-ion calculations a divergence-free perturbed field can matter.
With reg_flag=t the psi,theta grid output field is not divergence-free,
so a new output option reg_div_flag (default f) damps only the (n,m)
projection of (bwt,bwz) in gpout_xbrzphifun, which keeps div b = 0;
bwp, xi and all other outputs are unchanged.

gpout_arzphifun now takes the grad(zeta) components w(3,1:2) from the
chain rule through grad(psi) and grad(theta), as in gpeq.f, so that
curl A is consistent with the xbrzphifun field.

input/gpec.in template updated.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant