Skip to content

Merge upstream PrincetonUniversity/GPEC develop (54991242) into the fork - #10

Merged
krystophny merged 384 commits into
developfrom
sync/upstream-develop-20261102
Nov 2, 2026
Merged

krystophny merged 384 commits into
developfrom
sync/upstream-develop-20261102

Conversation

@krystophny

Copy link
Copy Markdown
Member

This PR records upstream ancestry with a real git merge --no-ff upstream/develop, with no rebase or squash. Upstream develop is at 5499124, release 1.6.0. #4 and #8 brought upstream's content in, but because they were squash-merged, git still counted 382 upstream commits as unmerged. After this merge, git rev-list HEAD..upstream/develop is 0.

What upstream brings in

No new content. Before the merge, git diff upstream/develop origin/develop showed only the fork's own additions, 8 files and +145 lines, and none of them removed anything:

Conflicts and resolutions

There was one conflict, in equil/direct.f, in the direct_run loop over itheta. It is left over from the squash. Upstream added that whole field-line sampling loop after the old merge base 1055d1b. The fork has the same loop plus #6's CALL direct_refine(ys(itheta,2),ys(itheta,0),psio*(1-sq%xs(ipsi))). Upstream never had that call and never removed it. I kept the fork side, and every line is within 72 columns.

make_test.yml and gpec/gpout.f merged automatically.

Content diff against current develop

git diff origin/develop 236ab2cd is empty. The merged tree is byte-identical to develop, so this PR only adds history.

Test evidence (local, gfortran 13, CI's make_gpec_gfortran FFLAGS, system NetCDF-Fortran)

  • cd install && make -j: all modules compiled successfully.
  • Soloviev ideal (dcon, gpec, stride, pentrc): pass
  • Soloviev kinetic (dcon, gpec, pentrc): pass
  • Soloviev resistive (rdcon, rmatch, gpec, stride): pass
  • Soloviev recon (dcon, gpec, check_recon.py): PASS, relative difference 6.97e-03 against a tolerance of 2e-02
  • pytest regression/test/test_direct_projected_radius.py regression/test/test_direct_error_message.py: 2 passed

Landing: the develop merge queue is squash-only, and squashing would drop the upstream ancestry again (as happened with #8), so this must land as a merge commit and is not enqueued.


Generated by Claude Code

jaesun57 and others added 30 commits July 20, 2026 14:14
(cherry picked from commit 5450935)
Add recon_flag3 path computing the |C|^2 decomposition of delta-W_p:
- gpeq_recon3: A=|Q|^2/mu0, B=2Re(V*.Q)/mu0 (V computed directly from
  the equilibrium current, not via C-Q), Cpar=K2 xi_n^2, Iperp term,
  plus optional per-theta surf_extra(0:mthsurf,22) field pack.
- gpout_recon3: 30-column surface output (gated by recon_out) and a
  C-free dW_p = 1/2(A+B+Iperp-dst1-dst3) via gpeq_dst (dst2 cancels).
- Unify recon1/2/3 terminal+log output (dW_p, dW_p-normalize, GPEC ep,
  DCON ep); fix recon1 DCON ep normalization to use dcon_fac.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 3f1f64b)
Add gpec_recon_terms_solX.out (per-point C-form energy densities:
c2_den, c2_psi/theta/zeta_den, dst1/2/3_den) from gpout_recon, the
recon1 analog of recon3_terms for 2 dW_p = C2/mu0 - dst1 - dst2 - dst3.

gpeq_epf / gpeq_dst gain optional per-theta density out-args; their
theta loops now run 0..mthsurf with an itheta<mthsurf integration
guard so the flux-surface integrals are unchanged.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
(cherry picked from commit f80b993)
main.tex section 4.1: add "Verification: C as an effective field" passage
(proof that (curl C).grad psi = 0; Eulerian vs Lagrangian frames; effective
current j_eff = curl C / mu0 by analogy with Ampere's law) and Figure 1
(C_eulerian_lagrangian_n1.png, force-added since *.png is gitignored). Also
reframe the V-component block as the correction term, drop section 6.1
(dcon_fac moved to the section-6 intro), and rename sections 4.1/4.2.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 73aaf92)
- A.11: replace the proportionality with an exact equality; expose the
  1/(B^2|grad psi|^2) prefactor and that the numerator vanishes by force
  balance mu0 j x B = mu0 p' grad psi.
- A.2 intro: correct which V-products enter H4. Only Y*.V, U*.V, and
  |V|^2 (with the -mu0 K scalar) form H4; X*.V feeds the C block and
  Z*.V = mu0 p' xi_psi'^* xi_psi feeds the E block.
- Remove redundant "Input Variables and Parameters" table.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 1ce5e3c)
Appendix A.1 presented V only through its components A_2,A_3, which read
backwards (the components are the *result* of expanding the physical field).
Fix the framing:

- State up front that V is the edge current rotated into the surface,
  V = (xi_psi/|grad psi|^2) mu0 j x grad psi -- a concrete (psi,theta)-grid
  field -- and that A_2,A_3 are just its Jacobian-weighted contravariant
  components.
- Add the scalar-triple-product rule (A.4):
  W . V = xi_psi (a mu0 j^zeta - b mu0 j^theta) for W = a v_2 + b v_3,
  i.e. the cofactor identity in coordinate-free form.
- Replace the repeated per-product cofactor derivations for X*.V, Z*.V, U*.V
  with one-line references to (A.4); drop the now-redundant cofactor
  intermediate equation.

No change to results; presentation only.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 69c6edc)
gpout_recon: dump gpec_recon_qvmn_sol*.out = covariant Q_i (bare delta-B = bv*_mn, captured right after gpeq_epf before gpeq_K/gpeq_dst). Effective field is cvmn = C_i = Q_i + V_i; cwmn is the contravariant J*C^i.

The quantitative QvsC figure built its (curl Q).grad(psi) panel from cwmn (contravariant) fed into the covariant curl operator -> spurious RMS 3.8e-2. Rebuilt from qvmn (covariant Q_i): RMS 2.1e-3 (Q pierces) vs 2.1e-12 (C tangent), matching the C-verify diagnostic; both panels now use covariant components as the caption states. Q-vs-C table + caption updated; figure swapped C_eulerian_lagrangian_n1.png -> QvsC_quant_n1.png (+ gitignore exception).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit fc716a4)
The xi_flag/xifile namelist entries + declarations were WIP and never wired up (idcon_inject_wt was never called). Removed from gpec.f, gpglobal.F, and input/gpec.in to avoid future conflicts.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 602eba2)
Bring 107 develop commits (RDCON/RMATCH resistive work, harvest removal,
positive-q sign conventions, install dependency regeneration) into the
recon energy-decomposition branch.

Conflict resolutions:
- gpglobal.F / gpec.f: keep both sides (recon covariant arrays, gpout_*
  types, recon_* flags + develop bwp_mn_rmatch, bpsi, use_res_spot,
  singthresh_slayer_inpr_prof)
- gpeq.f: keep recon's new subroutines (gpeq_c .. gpeq_recon3); take
  recon-side subprogram numbering
- out_ahg2msc / vac_memory: adopt develop's refactor (moved gpec -> equil
  module); dropped recon's stale gpec-side copies
- solovev gpec.in: take develop's example settings
- CLAUDE.md: keep deleted (recon removed it from VCS; local copy retained)

Builds clean: make gpec succeeds, gpec executable produced.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

# Conflicts:
#	gpec/gpec.f
#	gpec/gpeq.f
#	gpec/gpglobal.F
- remove unused gpout_shear/curvature/k bicube globals (allocated
  but never bicube_alloc'd, filled, read, or deallocated) and the
  banner text claiming bicube storage
- gpout_recon3 banner: B is computed independently from the
  equilibrium current, not derived as epf - A - Cpar - Iperp;
  identity_res is a genuine check (verified ~ -3e-10 on solovev)
- state that the Cpar = dst2 cross-check is not performed in-code
- move ep_index range checks to the top of recon2/3 (fail fast,
  matching gpout_recon); totals verified bit-identical
- document the psifac lower-bound handling difference between
  gpout_recon (clamp) and recon2/3 (skip)
- delete two dead declaration blocks (C-comparison diagnostics,
  unused *_cum accumulators)
- gpeq_c: comments said contravariant components are extracted by
  dividing by the Jacobian; no such division exists - the routine
  carries J Q^i / J C^i throughout (Jacobian-weighted convention)
- gpeq_cova header: metric%cs stores g_ij/J, so outputs are plain
  covariant, not Jacobian-weighted
- gpeq_cveri <-> gpout_recon inline copy cross-referenced
- fix wrong debug trace name in gpeq_c and gpeq_K exit message

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 3255c2f402f6afb8866e3d73aa9310099f8c0e56)
(cherry picked from commit 143bdaa1fd93b24dbdaf769eeea52711c92a73a4)
recon_flag3 drives the Bernstein |C|^2 decomposition (the recon3
diagnostic) but was missing from input/gpec.in; also note that
recon_out gates recon_flag2/recon_flag3 file output as well.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit cca31be7679547930599540020e542b47057bbcd)
(cherry picked from commit 47e86e5c7d5a2c8cab65973a8115fcd2cef9bee1)
Name the three GPEC_OUTPUT switches in the log-output section:
recon_out (file output gate), recon_int (spline/trapezoid psi
quadrature), cveri_flag (curl-C verify with PASS/WARN/SKIP against
1e-10), and label the cveri paragraph for the cross-reference.
Includes the current working revision of the derivation text.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 8ad60e848273cd9bbc870be8d721cad7cdd1c4ae)
(cherry picked from commit f89d399ed82155e708d54430f1f0670879a7295b)
…meter and set the wanted number of solution component then this will make Solution{sol_num}.bin file containing the {sol_num}th solution on Wt. If there aren't any sol_num then this will work just as same as before.

(cherry picked from commit 89e4516)
(cherry picked from commit 9294856)
…ializes sol_num to 1, allowing for better control of solution components.

(cherry picked from commit d975697)
(cherry picked from commit 3cbd047)
… index

(cherry picked from commit 25ddd6e)
(cherry picked from commit b5b3203)
Match input/match.in comment to match/match.f; sol_num=1 selects the most unstable DCON ideal eigenfunction. Addresses snu_offshoot review Q2.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
(cherry picked from commit 311922c)
(cherry picked from commit f998b3b)
- range-check sol_num against mpert (read from euler.bin at runtime)
  before indexing wt(:,sol_num); out-of-range values previously read
  past the array and wrote a garbage eigenfunction silently
- same check for mripple against mlow..mhigh on the ripple path
- remove the unconditional 'sol_num = ' debug write
- ripple_flag ignores sol_num, so write solutions.bin rather than a
  sol_num-labeled filename in that case

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
(cherry picked from commit 0bf08bf)
gpeq_sol/contra/cova/normal + gpeq_c were repeated inline at every
recon call site; bundle them as gpeq_prep_c(psi) and make the recon
kernels (gpeq_cveri, gpeq_epf, gpeq_dst, gpeq_epfterms) assume the
state is prepared. gpout_recon now calls gpeq_cveri directly instead
of carrying an inlined copy of the identity. Also fix a stray minus
in the cveri_flag comment of input/gpec.in.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
gpeq_shear, gpeq_curvature and gpeq_K assigned the gpglobal module
variables p1/q/q1 through host association. p1 in particular was
stored as p' while the rest of GPEC keeps mu0*p' there (idcon.f,
gpdiag.f). No reader consumes the stale values today, but shadow
them with routine-local declarations so the recon diagnostics leave
module state untouched. Results are bit-identical. Also note in
gpeq_c that the removed xno/bno and theta-space c2v* fields were
development-time verification aids.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
logan-nc and others added 17 commits October 24, 2026 19:11
…gauge-docs-20261019

DOCS - BUG FIX - State pmodb displacement gauges and prove consumer contracts
Keep the field-line diagnostic fix but defer the new source-extracting
test architecture until after the v1.6 release.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
gpout_xbnormal formed |grad psi| from the shared r(:) without setting
it, so each surface used whatever surface the previous caller left
there: usually one surface behind, and on the first surface whatever
gpout_singfld left. Set r(itheta) from the surface's own rs, as the
gpeq_* routines do. xno, bno and bwp move by up to 0.4% near the edge
and no longer depend on singfld_flag.

Supersedes PrincetonUniversity#353.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE
…n the estimate, and stop on bad layer inputs

- sglobal.f: lnLamb was an uninitialized variable set only in params(),
  so gpec_slayer (which never calls params) ran with lnLamb = 0 and
  produced infinite S and Q. Initialize it to 17, the old PARAMETER
  value.
- slayer.f: the growth-rate estimate loop set only D_norm per surface;
  w_der_del_s also reads tau and P_tor, which still held the last
  surface's values. Set them per surface as the matching loop does. On
  the DIII-D example the 2/1-5/1 estimates change by 12-25%.
- params.f: stop with a message naming the bad value and the surface
  when n_e, t_e or t_i is not positive, when the shear is zero, or
  when chi_perp is not positive where the critical Delta-prime needs
  it. These previously gave NaN or Inf, or an unexplained
  non-convergence stop.
- TARGETS.inc: match only v* tags in the slayer version rule, as the
  other codes do. On a full clone the tps-pin-* tags gave a version
  string past column 72, and slayer_netcdf.f failed to compile.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE
…ugfix/xbnormal-surface-r

GPEC - BUG FIX - Use each surface's own R for |grad psi| in xbnormal

This counts as reviewed as it is itself the review product of 353
…ine-diagnostic-upstream-20261023

BUG FIX - Preserve complete field-line integration failure diagnostics
…elist mode, use per-surface shear, ds and Q in the Br threshold, drop dc_type lar, and carry |s| and n into every path

Closes the six items left open in the review of PrincetonUniversity#173.

- slayer.f: the namelist has no silent defaults since this branch
  (zero means unset, or "compute" for inQ, inQ_e, inQ_i). A run that
  needs an unset input now stops with the variable, the value read and
  its v1.5 default, instead of running on zeros (inpr=0 gave LSODE
  atol=0, inds=0 gave 0/0 in the Riccati boundary condition). The
  dimensional inputs are checked when params() runs on the namelist,
  the Cole inputs when the single-surface riccati() block runs. That
  block is now also skipped for br_th_flag, which computes its own
  deltas per surface like the growth-rate modes.
- slayer.f: gamma_fac was never set in namelist mode, so the estimate
  was 0. It is now rs*Re(delta_prime)/tau_r after params(), as
  build_inputs forms it from the r_s-referenced Delta_prime. The
  template with delta_prime=(5,0) estimates 102.7 Hz.
- sglobal.f, growthrates.f, layerinputs.f, slayer.f: slayer_inputs_type
  carries shear_arr, ds_arr and Qexb_arr. The Br threshold loop used
  the namelist sval on every surface, Q0 = Q_e where the torque balance
  needs the ExB Q, and the Fitzpatrick D_norm in the Cole ds slot of
  riccati(). It now uses each surface's shear, ds and Q, as the other
  riccati() call sites do. On the DIII-D example with read_eq=t it
  gives finite thresholds on all five surfaces.
- params.f, slayer.f: dc_type='lar' is removed. It was the rfitzp
  formula without the r_s normalization (lar = rfitzp/r_s), so it was
  dimensionally inconsistent and, corrected, a duplicate. It stops with
  a message like 'toroidal'.
- gslayer.f: the GPEC path uses ABS(sval) in the Alfven time and n in
  Q_e, Q_i (GPEC passes the n = 1 diamagnetic frequencies), the same
  fixes ffa95f1 and d43cfb5 made in params.f and build_inputs.
  params.f: the namelist omega_e, omega_i carry n too. No change at
  n = 1 with positive shear.
- DIIID README: the coupled scan needs msing_max=2 or a wider
  scan_width, since the Doppler shifts put the 2/1 outside the shipped
  window. With msing_max=2 the coupled root is omega = 30.93 krad/s,
  gamma = -379.6 /s.

The DIII-D uncoupled run, a10 and the develop template give identical
screen output and NetCDF before and after.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…nput report, and say plainly where GPEC's n enters

Follow-ups from an independent review of 2100323:

- qval is needed with l_n and l_t for the namelist omega_e, omega_i;
  without it params() stopped in slayer_d_norm with a Q_e = Q_i message
  instead of naming qval. It is now in the report (v1.5 default 2.0).
- The IF (params_flag) block calls params() on the namelist values
  whatever read_eq is, and params_flag defaults to true, so the
  dimensional inputs are now checked whenever params_flag is set, not
  only when read_eq is false. read_eq=t with params_flag=t and no
  density now reports n_e instead of failing inside params().
- The report's hint about read_eq applies to the dimensional inputs
  only; reworded. The Cole-input comment no longer says the legacy
  scans are skipped by est_gamma_flag (they are not).
- gslayer.f, params.f comments: GPEC hands gpec_slayer the bare
  diamagnetic frequencies with no toroidal mode number attached, and
  GPEC's own nn multiplies them there; slayer.in is not read by GPEC.

The DIII-D run, a10 and the develop template are unchanged.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…an, guard the Br threshold against NaN, and report Q_sol at the right grid point

Closes the remaining points from the independent review of PrincetonUniversity#173 so
that 1.6.0 leaves no SLAYER follow-ups behind.

- gslayer.f, params.f: the normalised ExB Q now carries the toroidal
  mode number like Q_e and Q_i already did (Q = Qconv*n*omega_E), the
  Doppler frequency the layer sees and the convention of the coupled
  determinant's shift. Before, the torque balance 2*pr*(Q0 - Q)
  compared an n-free Q0 with n-scaled Q_e for n >= 2. gpec_slayer's
  omega_sol is divided by n again so it stays a fluid rotation
  frequency. GPEC passes the bare frequencies and its own nn;
  slayer.in is never read by GPEC. No change at n = 1.
- slayer.f: the stability scan with Pperp_Ptor_flag=f passed the
  Fitzpatrick D_norm in riccati()'s Cole ds slot; it now gets the
  surface's ds, as the Br threshold does since 2100323.
- slayer.f: the single-surface baseline riccati() is also skipped for
  input_flag, whose per-surface results never used it, and the Cole
  inputs are required whenever a legacy scan will run.
- slayer.f: the Br threshold (bal_flag and br_th_flag) ignores
  non-finite torque-balance entries and reports 0 with a warning when
  the balance has no positive maximum, as gpec_slayer does, instead of
  NaN. The a10 example now prints br_th = 0 with the warning.
- slayer.f: MAXLOC counts from 1 while the scan arrays start at 0, so
  every reported Q_sol was one grid step above the maximum. Fixed at
  all seven sites. gpec_slayer uses the same point only as the centre
  of its refinement bracket, so it is left as is and GPEC's output is
  unchanged.
- delta.f: riccati()'s inds argument is ds, the normalised ion sound
  radius, not a shear.

The DIII-D run and the develop template are unchanged; the a10
example's Q_sol moves from 0.05 to 0.0 and its br_th from NaN to 0.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…layer_growthrate

Growth rate calculation functionality added to gslayer.f
…ugfix/pentrc-wmat-conjugation

PENTRC - BUG FIX - Conjugate the bounce-averaged W vectors so the kinetic EL matrices see the physical perturbation
…elease-1.6.0

Merge release-1.6.0 into develop
Records upstream ancestry lost when #4 and #8 were squash-merged. The
content of upstream develop 5499124 is already in the fork, so the
merged tree equals the previous fork develop.

Conflict: equil/direct.f (direct_run loop). Upstream added the
field-line sampling loop after the merge base; the fork carries the same
loop plus the direct_refine call from #6. Kept the fork side.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01H2kuiYoFvpcfgPswUWX5Rv
@krystophny
krystophny added this pull request to the merge queue Nov 2, 2026
@github-merge-queue
github-merge-queue Bot removed this pull request from the merge queue due to failed status checks Nov 2, 2026
@krystophny
krystophny added this pull request to the merge queue Nov 2, 2026
@github-merge-queue
github-merge-queue Bot removed this pull request from the merge queue due to failed status checks Nov 2, 2026
@krystophny
krystophny added this pull request to the merge queue Nov 2, 2026
@github-merge-queue
github-merge-queue Bot removed this pull request from the merge queue due to failed status checks Nov 2, 2026
krystophny added a commit to itpplasma/.github that referenced this pull request Nov 2, 2026
A fork that merges its parent repository inherits commits whose
identities, file sizes and contents it cannot change. In a GitHub fork
the workflow now fetches the parent's branches into refs/upstream/* and
the script excludes commits reachable from them (--upstream-ref).
Identities and secrets are checked on the fork's own commits only. A
(path, blob) pair introduced by a parent commit stays exempt where it
reappears, for example in the fork's merge commit; any fork change to
such a path, and any new path, is checked as before.

On itpplasma/GPEC#10 (develop merged with PrincetonUniversity/GPEC
develop) the gate previously failed on 57.6 MB of parent blobs and three
parent identities; with the parent refs only the merge commit is
checked, and it passes. The job timeout rises to 10 min for the parent
fetch.


Claude-Session: https://claude.ai/code/session_01H2kuiWu67fzfh8w4Fkby3P

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
@krystophny
krystophny added this pull request to the merge queue Nov 2, 2026
Merged via the queue into develop with commit 2d4443b Nov 2, 2026
6 checks passed
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.

9 participants