[HengyuLi-Ozaki-lab]

The rate-path PR asked for in the triage on #236: one commit per defect, issue number in the subject, on develop 4b5463e4. The fixes change trx/libnf.f90 and trx/trpnf.f90; of the three commits added after your review, one also touches trx/trmenu.f90 (a printed unit) and one CHANGELOG.md.

commit	item	change
39dfc77c	#234-2	pm_local is REAL(rkind) (it truncated a mass of order 1e-27 kg to 0)
26cd2f3e	#234	module-level IMPLICIT NONE in libnf_local and libnf; compiles cleanly (#234-3 does not apply to develop, as you said)
c0dd2721	#236	sigmav_nf = sigmav_nf*1.D-6 at the single computed exit
31e65493	#235-1 (#229-3)	PA, PZ, COULOG of the slowing-down loop take nsp, not the stale ns
ef975281	#235-2 (#229-2)	tr_pnf zeroes PNF_NSNNFNR with its siblings
293899fe	#229-1	tr_pnf zeroes PNFIN_NSNNFNR too (whether it should carry a value is left open)
3953ef52	#235-3	CASE(1,12) / CASE(2,3) / CASE(4,14), and libnb -> libnf in the message (14 is refused again by 428c0307)
1cae5566	#235-3	nsp_idnf(id_nf_dd3): He4 for the reduced variants 12 and 14, He3 for 2, 3 and 4 (its own commit, so it can be dropped)
5ab61000	#229-4	the dd2 block sets enn_idnf(id_nf_dd2), not dd3's
5c7b87d4	#268 (new)	VF = SQRT(2*eng_nnf/(PA(nsp)*AMP)): eng_nnf is already in joules (eng_idnf(id_nf_dt)=3.5D3*RKEV), and the *RKEV on this line applied it twice
428c0307	#235-3 (review)	model_pnf=14 gets its own CASE(14) and stops with "no He5 species can be defined"; CASE(4) is as on develop
c4bca92b	#236 (review)	sigmav_nf's unit, m^3/s, in its two comments and in the menu's B output
8daf0d6c	(review)	a CHANGELOG.md Fixed entry (for this PR only since b9ae2485)
b9ae2485	(second review)	the entry covers this PR only (its #270 sentences moved to #270's d2bbfcca); model_pnf=14 marked not available in the model tables of trx/libnf.f90 and trx/trinit.f90; the CASE(14) comment reworded

model_pnf=14 stays refused (428c0307, from your review). I had written in this description and on #236 that He3 cannot be defined from the menu, and concluded that 14 only reaches its presence check. That was wrong: trmain reads plparm, whose &PL takes NPA and PA. With H and He3 defined there (NSMAX=6), CASE(4,14) let 14 through on NS_He5 = 7, and tr_pnf stopped at trpnf.f90:88 with "Index '7' of dimension 1 of array 'snf_nsnnfnr' above upper bound of 6" (reproduced here). 14 now has its own CASE(14) and stops with "XX Error: libnf: model_pnf=14 is not available: no He5 species can be defined". CASE(4) is as on develop: its same He5 hole goes to the lead, with the species set 14 is meant to have and its remaps, which do not agree with each other (the32 treats it as reduced; dd2, dhe32, the35, the36 as the full variants).

The dd3 remap. The header gives 12 and 14 the same DD3 line, D + D -> He3 + n ! He4 + n, so both collapse the He3 product onto He4. On develop the remap gives He3 to 12 and He4 to every other variant; 1cae5566 gives He4 to 12 and 14 and He3 to 2, 3 and 4. So 12 changes from He3 to He4 and 2, 3, 4 from He4 to He3: for those three a choice that changes results, as you note, and the reason the commit stands alone. With 14 refused, its branch there is not reached.

5c7b87d4 was in no report; it is now #268. It is the fifth correction on our reference branch, and it sits on the same line as #235-1. Without it VF is 0.16 m/s instead of 1.3e7 m/s (a factor sqrt(RKEV)), so VF/VCR is about 4e-8 to 1.3e-7 on this deck instead of about 3 to 10. There 1-HY(VF/VCR) in TAUF should be of order 1e-22, but in double precision it is round-off of order 0.01 to 0.1 (and -0.69 at 1e-8), so TAUF is meaningless and can be negative. The alpha tail energy WTAILT of the intermediate rows below therefore means nothing; with the commit it is 4.27e-2 MJ. It is a commit of its own, so that it can be dropped if you prefer.

One effect on #226. That branch still declares PUBLIC sigmav_nfb_int (trx/libnf.f90:117) with nothing behind it. Under 26cd2f3e's module-level IMPLICIT NONE the line no longer compiles. Git shows no conflict there and CI does not build trx, so nothing would flag it; I will drop the line when #226 merges develop.

Left out, as asked: #234-1 (the Duane columns), #238's secondary items, and the STOP -> ierr conversion (#142).

How it was built and run

develop's CI does not build trx, so this is how it was run here: macOS arm64, gfortran 15 (MacPorts), the make.header flags with run-time checks (-fbounds-check -fcheck=all -ffpe-trap=invalid,zero,overflow), GSAF 4.16 and X11 for the graphics that tr2 links, and mtxp on MUMPS with MPICH; make -C trx; trx/tr2 < deck in an empty directory. The deck is #238's 10 keV D-T case (MODELG=2, e/D/T/He4, model_pnf=1), with w 1 x added before q so that TRPRNT type 1 prints the global results. The same deck was also run with model_pnf=0, 12 and 14, and a 1 keV variant of it. The decks draw nothing, so their GSAF files are empty and are not compared. The variants that need He3 (2, 3, 4) need it defined in a plparm file (&TR takes no NPA). With H and He3 there and NSMAX=6, model_pnf 0, 1, 2, 3 and 12 all stop in the first step at trexec.f90:533, "Index '0' of dimension 2 of array 'b' below lower bound of 1", on develop and on this branch alike, so no run with He3 completes here (2 and 3 now pass their check and dd3 gives them He3, but there is no result for them); now #275. model_pnf=14 was run both ways: on this deck, and with that plparm (below).

The 10 keV deck

0
f
trx_fus_dt_hot.gs
c
MODELG=2
RR=8.481D0
RA=2.574D0
RKAP=1.816D0
RDLT=0.3478D0
BB=5.953D0
NSMAX=4
PN=0.1D0,0.045D0,0.045D0,0.005D0
PNS=0.01D0,0.0045D0,0.0045D0,0.0005D0
PT=10.0D0,10.0D0,10.0D0,10.0D0
PTS=1.0D0,1.0D0,1.0D0,1.0D0
PROFN2=0.15D0,0.15D0,0.15D0,0.15D0
model_pnf=1
DT=0.02D0
NRMAX=50
NTMAX=5
NTSTEP=5
RIPS=2.D0
RIPE=3.D0
r
w
1
x
q

PROFN2 is per species in trx, so all four values are written out.

What each commit does to the runs

trx rebuilt from clean at develop and at each commit (WPT and WTAILT in MJ; where a commit was later rewritten in its message or CHANGELOG.md only, the runs are those of the same trx sources):

at	10 keV, model_pnf=1	model_pnf=12	model_pnf=14
develop	aborts with negative temperatures at nr=31 (reported time 0.040 s; eight diagnostics, XX TRXTOA3 and XX TRXTOA4)	stops at trpnf.f90:86: "Index '0' of dimension 1 of array 'snf_nsnnfnr' below lower bound of 1" (the dd3 product is NS_He3, which is 0)	XX Error libnb: undefined model_pnf, STOP
39dfc77c, 26cd2f3e	the same	the same	the same
c0dd2721 (#236)	runs: WPT 28.71, WTAILT 8.28e-2	the same	the same
31e65493 (#235-1)	WPT 31.07, WTAILT 2.19 (the fusion power still accumulates)	the same	the same
ef975281 (#235-2)	WPT 28.62, WTAILT -3.19e-3	the same	the same
293899fe (#229-1)	unchanged	the same	the same
3953ef52 (#235-3)	unchanged	the same	XX Error: libnf: NS_H=0 (the check is reached; the deck has no H)
1cae5566 (dd3)	unchanged	runs: WPT 28.70; dd2 and dd3 give He4, dd1 T	the same
5ab61000 (#229-4)	unchanged	unchanged	the same
5c7b87d4 (VF)	WPT 28.65, WTAILT 4.27e-2	WPT 28.65, WTAILT 4.30e-2	the same
428c0307 (14 refused)	unchanged	unchanged	XX Error: libnf: model_pnf=14 is not available: no He5 species can be defined
c4bca92b, 8daf0d6c, b9ae2485 (comments, the B print and CHANGELOG.md only)	unchanged	unchanged	the same

model_pnf=14 with H and He3 defined in plparm (NSMAX=6): on develop XX Error libnb: undefined model_pnf; at 5c7b87d4 the stop at trpnf.f90:88 quoted above; from 428c0307 on the new message.

The plparm file and the deck for it

plparm, in the run directory:

 &PL
  NSMAX=6
  NPA(5)=1, PA(5)=1.D0, PZ(5)=1.D0
  NPA(6)=2, PA(6)=3.D0, PZ(6)=2.D0
 /

The deck is the 10 keV deck above with these lines in place of its NSMAX to model_pnf lines:

NSMAX=6
PN=0.1D0,0.044D0,0.044D0,0.005D0,0.001D0,0.0005D0
PNS=0.01D0,0.0044D0,0.0044D0,0.0005D0,0.0001D0,0.00005D0
PT=10.0D0,10.0D0,10.0D0,10.0D0,10.0D0,10.0D0
PTS=1.0D0,1.0D0,1.0D0,1.0D0,1.0D0,1.0D0
PROFN2=0.15D0,0.15D0,0.15D0,0.15D0,0.15D0,0.15D0
model_pnf=14

The B command now prints sigmav [m^3/s]= 1.1000000000000000E-022 for D-T at 10 keV.

With model_pnf=0 every commit prints the same output as develop apart from the CPU and elapsed-time lines. On the 1 keV deck TAUE1 goes from 3.16e-7 s on develop to 0.267 s at c0dd2721 and 1.164 s from ef975281 on. Two consecutive runs of each deck give the same output apart from those time lines.

PINT, TAUE1 and QF are still corrupted by #238 on this branch: on the 10 keV deck PINT is 1.109e5 where it should be about 0.54 MW, so PINT and TAUE1 are off by about 2e5 and QF by 1e6. That is the separate PR. A second R in one process ends trx with "Attempt to DEALLOCATE unallocated 'gpnb'" (trx/trcomm.f90:996), on develop as well, so there is no in-process two-run check here; not touched in this PR. No CI test is added, since CI does not build trx and tr2 needs the graphics library.

🤖 Generated with Claude Code

[k-yoshimi]

Review through the full cycle (Codex + in-house domain, test-oracle and verifier lenses + docs/workflow review), with trx/tr2 built here at develop, this head, #270 and the two merged (gfortran 15.2, the make.header run-time checks, graphics stubbed with tr/tr_static_stubs.f90).

Reproduction. Your deck with model_pnf=0 gives output identical to develop on this head, on #270 and on both merged (CPU-time lines stripped), and identical across repeated runs. With model_pnf=1: WPT 28.65 MJ, WTAILT 4.265e-2 MJ on this head; merged with #270, PINT 0.5393 MW, QF 1.294, TAUE1 53.13 s, TAUE2 1.518 s — your numbers to four digits. Every fix was traced independently: the tables are cm³/s (the trm sibling and trx/trsigmavnf.f90 apply the same 1e-6), eng_nnf is in J, the stale ns is NSMAX+1, tr_pnf runs at least three times per step so the accumulation was real, IMPLICIT NONE compiles clean, model_pnf=13 was already unreachable.

One change requested — 3953ef52, the CASE(4,14) label. Before this PR model_pnf=14 reached CASE DEFAULT and stopped. After it, 14 passes the presence check whenever D, T, He4, H and He3 are present, because NS_He5 is never a discovered index: pl/plinit.f90:122 sets it to the literal 7, tr_prep_ns (trx/trprep.f90:122-160) neither resets nor matches it, and pl_prep_ns (pl/plprep.f90:50-85) refuses an NPA=2, PA=5 species outright. So for 14 nsp_idnf(id_nf_the35) is 7 (libnf.f90:311) and tr_pnf writes SNF_NSNNFNR(7,…) / PNF_NSNNFNR(7,…) past the (NSMAX,NNFMAX,NRMAX) allocation (trx/trcomm.f90:686,688), and trexec.f90:568/643/720 read it back. Reachable from tr2 with a plparm file that defines He3 (NPA=2, PA=3), i.e. exactly the composition 14 is documented for. Please keep 14 refused for now: drop the ,14 (so it stops at CASE DEFAULT as on develop), or give it its own CASE(14) with a message saying He5 cannot be defined. The dd3 commit 1cae5566 can stay; its 14 branch is then unreachable. Note for the lead, kept separate from the mechanical fixes: that commit also moves the DD3 product for model_pnf=2/3/4 from He4 to He3, which is the species the header documents and which the domain reviewer endorsed, but it is a result-changing physics choice of the same class as the Duane mapping; it is its own commit so it can be dropped if the lead prefers. Which species set 14 is meant to have (the header says nsmax=4, trinit.f90:490 says 6, the check needs 7) and whether model_pnf=4 — which has the same He5 hole today — should be refused too, goes to the lead with the other items you listed.

Found on the way, pre-existing, filed as #271 and #272 rather than asked of this PR: the natural cubic spline of the linear svnf_* values against log10(T) returns negative reactivities for D-He3 (2–4.5 keV) and T-He3, and a D-T value 5× too high at 1.5 keV and 11× too low at 2.5 keV (I reproduced this with lib/libspl1d.f90 and the tables; trm/libsigma.f90 splines log10(σv) instead); and the two-entry reactions (DD1/DD2, DHe31/32, THe3x) each subtract SNF from both reactant slots, so D is depleted once per entry rather than once per reaction. Both change model_pnf>0 results again once fixed; your 1e-10 oracle will be useful there.

LOW, fold in if you like: the sigmav_nf header comment could state the return unit (m³/s), and trmenu.f90:157 prints sigma= with no unit, now in m³/s. A CHANGELOG.md Fixed entry covering this PR and #270 (results of model_pnf>0 runs change) would be welcome; trx is outside the file's stated Phase L scope, so your call.

After the one change: approve from my side; merge needs the lead's confirmation as well, per the plan on #236.

[HengyuLi-Ozaki-lab]

Thank you. The change you asked for is done, and so are both LOWs: three new commits on top of 5c7b87d4, nothing before it rewritten.

• 428c0307: model_pnf=14 has its own CASE(14) and stops with XX Error: libnf: model_pnf=14 is not available: no He5 species can be defined; CASE(4) is as on develop.
• c4bca92b: sigmav_nf's PUBLIC line and header comment give m^3/s, and the B command prints sigmav [m^3/s]=.
• 8daf0d6c: one CHANGELOG.md Fixed entry for this PR and fix(trx): the ten per-channel power integrals of TRGLOB convert W to MW (#238) #270. fix(trx): the ten per-channel power integrals of TRGLOB convert W to MW (#238) #270 itself is unchanged. If the lead drops the dd3 commit or 5c7b87d4 (trx/trpnf.f90: the fusion-product speed VF applies RKEV twice (eng_nnf is already in joules) #268), or fix(trx): the ten per-channel power integrals of TRGLOB convert W to MW (#238) #270 does not go in, the entry's words on the D-D branch to He3, on RKEV, or on TRGLOB and the NBI powers go with it. On NBI and RF it says what I measured rather than a 1e6 for PINT and TAUE*: with model_nnb=1, PNBIN=10 and with PECIN=5 on trx/trrslt.f90: ten per-channel power volume integrals omit the W->MW conversion their legacy siblings have; PNF_TOT/PNB_TOT/PRF_TOT enter PINT 1e6 too large #238's deck with fusion off, develop and fix(trx): the ten per-channel power integrals of TRGLOB convert W to MW (#238) #270 give the same PINT and TAUE1, and PNB_TOT and PRF_TOT are 0 on both, because the arrays they integrate are never filled (trx/trrslt.f90: ten per-channel power volume integrals omit the W->MW conversion their legacy siblings have; PNF_TOT/PNB_TOT/PRF_TOT enter PINT 1e6 too large #238); only the NBI run's printed PBCLE to PBCLA change, each divided by 1e6 (PBCLE from 5.123e5 to 0.5123); in the EC run they are 0 on both.

You were right about 14, and my description and my comment on #236 were wrong: He3 can be defined, through the plparm file that trmain reads. I reproduced your case with H and He3 in plparm and NSMAX=6: at 5c7b87d4, tr_pnf stops at trpnf.f90:88 with "Index '7' of dimension 1 of array 'snf_nsnnfnr' above upper bound of 6"; from 428c0307 on, the run stops with the new message. I have corrected the description: the 14 paragraph, the He3 sentence, and the rows for the new commits, with the plparm file used. With that plparm, every model_pnf I tried (0, 1, 2, 3, 12) stops in the first step at trexec.f90:533 on develop as on this branch, so no run with He3 completes; the description notes it for the lead.

Runs at the new head, trx rebuilt from clean: the 10 keV deck, the 1 keV deck and model_pnf=12 print the same output as at 5c7b87d4, and with model_pnf=0 the same as develop, in both cases apart from the CPU and elapsed-time lines. libnf.f90 compiles with no warning under -Wall -Wextra -Wconversion. The branch merges without conflict into develop ea9806ac and with #270.

The lead's items stay as you listed them: the species set 14 is meant to have, model_pnf=4's same He5 hole, and the dd3 choice for 2, 3 and 4. When #271 and #272 are worked on, our reference branch can give numbers for them.

[k-yoshimi]

Second pass over the delta 5c7b87d4..8daf0d6c (Codex + in-house domain and verifier lenses + docs/workflow review), with trx/tr2 rebuilt from this head the same way as before (gfortran 15.2, bounds and FPE checks, graphics stubbed).

The code change is accepted. The species check's SELECT CASE is well formed ((1,12), (2,3), (4), (14), DEFAULT), CASE(4) is as on develop, and the comment in CASE(14) is right on every point we checked (pl/plinit.f90:122, trx/trprep.f90:114-164, pl/plprep.f90). The unit annotations match what sigmav_nf returns (one *1.D-6 at trx/libnf.f90:443, tables in cm^3/s). Nothing parses the old sigma= label.

Reproduction. The 10 keV deck, the 1 keV deck and the fusion-off deck print the same output as at 5c7b87d4 apart from the CPU and elapsed-time lines. With your plparm (H and He3, NSMAX=6) and model_pnf=14 the run stops with XX Error: libnf: model_pnf=14 is not available: no He5 species can be defined, as you describe. With the same plparm and model_pnf=1 the run stops at trexec.f90:533 (Index '0' of dimension 2 of array 'b', NEA(NS,2)=0 for the species without an equation) on develop's tr2 as well, so that is pre-existing; It is now #275, so the lead's item has a number.

One thing to change, text only — CHANGELOG.md. All three reviewers raised it independently and the verifier confirmed it against the code: the entry is under Fixed on this branch, but two of its sentences describe trx/trrslt.f90, which this branch does not touch. In this checkout the ten TRGLOB integrals (trx/trrslt.f90:167-328) still have no /1.D6, and PBCLE..PBCLA (PNBCL_NS, printed at :930) are still in W. If #269 merges before #270, or #270 is held, develop's changelog says fixes exist that do not. Please:

• move the sentence "Ten power integrals of TRGLOB stayed in W …" and the sentence "In NBI runs the printed beam collisional powers (PBCLE to PBCLA) are now in MW; PNB_TOT and PRF_TOT …" into an entry on #270 (a docs-only commit there; we will re-check it without a rebuild), and tag this entry (#269; …);
• "never set PNFIN_NSNNFNR": the only write on this branch is the zeroing at trx/trpnf.f90:70, so PNFIN_TOT is still always 0. "is now zeroed each call" says what the code does.

Merging #270 first (or together) would also make the text true, and if the maintainers prefer that we can take it, but the split is cleaner whatever the lead decides.

LOW, same commit if you like.

• The parameter tables in trx/libnf.f90:37-44 and trx/trinit.f90:490-497 still list model_pnf=14 as a model; one word there ("not available: no He5 species") saves a user a run.
• the35's product in the CASE(14) comment reads as a typo for "the THe35 product"; and "past the arrays when NSMAX < 7" is true of the (NSMAX,NNFMAX,NRMAX) arrays (SNF_/PNF_/PNFIN_/PNFCL_NSNNFNR), not of RN/RT/PA/PZ, which are sized NSTM/NSM.

Recorded, not requested. model_pnf=4 has the same NS_He5 hole (the guard at :202 can never fire on NS_He5=7); that, the species set 14 is meant to have, and the dd3 commit stay on the lead's list as before. PNFCL_NSNNFNR is likewise only zeroed, so PFCLE..PFCLA are always 0 (same class as #238). The refusal is a STOP like its neighbours; trx/libnf.f90 is linked into tr2 only, not into any shared library, so #142 does not apply.

Verdict: request changes for the changelog text only; approve follows the push, no rebuild needed on our side for a docs commit.

[HengyuLi-Ozaki-lab]

Done in b9ae2485, docs only (no code change, so no rebuild needed):

• CHANGELOG: this entry is now tagged (#269; issues #229, #234, #235, #236, #268) and describes this PR only.
  • The sentences on the ten TRGLOB integrals and on the NBI powers moved to an entry on fix(trx): the ten per-channel power integrals of TRGLOB convert W to MW (#238) #270 (d2bbfcca). It sits after the second [Unreleased] Fixed entry, so fix(trx): the fusion rate path, one commit per defect (#234, #235, #236, #229, #268) #269, fix(trx): the ten per-channel power integrals of TRGLOB convert W to MW (#238) #270, feat(fp_mcp): expose set_param_str tool for KNAMEQ / MODELG=3 equilibrium load #231, fix(eq): eq_save refuses after a set of tot's eq: parameters too #267 and Reconcile bpsi/develop base into develop (bpsi↔fork merge baseline, CI-green) #226 do not conflict over it (checked with git merge-tree).
  • "never set PNFIN_NSNNFNR" now reads what the code does: PNF_NSNNFNR and PNFIN_NSNNFNR are zeroed on each call, and nothing else writes PNFIN_NSNNFNR, so its totals stay 0.
• The two LOWs: the model tables in trx/libnf.f90 and trx/trinit.f90 mark model_pnf=14 "not available: no He5 species can be defined". The CASE(14) comment names the THe35 product, and tr_pnf's (NSMAX,NNFMAX,NRMAX) arrays as the ones it would overrun.

libnf.f90 still compiles with no warning under -Wall -Wextra -Wconversion. The branch merges without conflict into develop and with #270, #231 and #267; with #226 only trx/libnf.f90 conflicts, as before. Thank you for filing #275.

[k-yoshimi]

b9ae2485 checked: the diff against 8daf0d6c touches CHANGELOG.md and comment lines only (no non-comment line of trx/libnf.f90 or trx/trinit.f90 changes), so the tr2 built from 8daf0d6c and its runs stand. The entry now describes this branch: every sentence in it is something trx/libnf.f90 or trx/trpnf.f90 does here, the PNFIN_NSNNFNR sentence says what the code does, and the two TRGLOB sentences are on #270 (d2bbfcca, approved there). The model tables mark 14, the CASE(14) comment names the THe35 product and the (NSMAX,NNFMAX,NRMAX) arrays. CI 9/9 on this head. git merge-tree here agrees with you: clean against develop, #270, #231 and #267; trx/libnf.f90 conflicts with #226, which is why #269 should merge before #226.

Verdict: approve. Merge together with #270 (either order now that each carries its own changelog entry), on the lead's confirmation per the plan on #236; the lead's items (species set of 14 and 4, NS_He5, the dd3 choice, #275) are unchanged.