diff --git a/CHANGELOG.md b/CHANGELOG.md index b27ede66a..40736285c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -96,6 +96,20 @@ The format is loosely based on [Keep a Changelog](https://keepachangelog.com/). ### Fixed +- **trx: the fusion reaction path** (#269; issues #229, #234, #235, #236, + #268). `sigmav_nf` returned its tables in cm^3/s where `tr_pnf` takes + m^3/s, so every `model_pnf` rate was 1e6 too large (#238's 10 keV D-T + deck stopped with negative temperatures). `tr_pnf` took the slowing-down + species from a stale loop index, applied `RKEV` twice to the product's + speed (so `TAUF` was round-off) and added `PNF_NSNNFNR` up over its + calls; `PNF_NSNNFNR` and `PNFIN_NSNNFNR` are now zeroed on each call + (nothing else writes `PNFIN_NSNNFNR`, so its totals stay 0). Also in + `libnf`: `pm_local` is a `REAL`, both modules are `IMPLICIT NONE`, the + D-D branch to He3 gives He3 for `model_pnf` 2, 3 and 4 and He4 for the + reduced 12 (which ran past an array bound before), the dd2 block sets + its own `enn_idnf`, and the species check refuses `model_pnf=14` with + its own message, since no He5 species can be defined. Results of every + `model_pnf>0` run change. - **MCP servers: the `mcp` dependency is pinned to `>=1.2,<2`** in the seven server `pyproject.toml` files, the CI workflow, the servers' own "not installed" message, the READMEs, the Sphinx MCP pages, the diff --git a/trx/libnf.f90 b/trx/libnf.f90 index b0dad6247..314ee7af7 100644 --- a/trx/libnf.f90 +++ b/trx/libnf.f90 @@ -2,13 +2,15 @@ MODULE libnf_local USE bpsd_kinds - INTEGER:: id_nf_local,pm_local - REAL(rkind):: temperature_local + IMPLICIT NONE + INTEGER:: id_nf_local + REAL(rkind):: pm_local,temperature_local END MODULE libnf_local MODULE libnf USE bpsd_kinds USE bpsd_constants + IMPLICIT NONE ! Fusion model ! model_pnf=0 : no fusion reaction @@ -40,6 +42,7 @@ MODULE libnf ! T + He3 -> He4 + p + n THe31 THe32 ! T + He3 -> He4 + D THe33 THe34 ! T + He3 -> He5 + p ! He4 + p THe35 THe36 + ! not available: no He5 species can be defined ! Fusion reaction id INTEGER,PARAMETER,PUBLIC:: id_nf_DT= 1 ! D + T -> + n @@ -112,7 +115,7 @@ MODULE libnf PUBLIC set_usigmav_nf ! set_usigmav_nf PUBLIC sigma_nf ! sigma_nf(id_nf,energy) - PUBLIC sigmav_nf ! sigmav_nf(id_nf,temperature) + PUBLIC sigmav_nf ! sigmav_nf(id_nf,temperature) [m^3/s] CONTAINS @@ -180,7 +183,7 @@ SUBROUTINE set_usigmav_nf END SELECT SELECT CASE(model_pnf) - CASE(1,12,13) + CASE(1,12) IF(NS_D*NS_T*NS_He4.EQ.0) THEN IF(NS_D.EQ.0) WRITE(6,*) 'XX Error: libnf: NS_D=0' IF(NS_T.EQ.0) WRITE(6,*) 'XX Error: libnf: NS_T=0' @@ -206,8 +209,16 @@ SUBROUTINE set_usigmav_nf IF(NS_He5.EQ.0) WRITE(6,*) 'XX Error: libnf: NS_He5=0' STOP END IF + CASE(14) + ! NS_He5 is not a species found in the run: pl_init sets it to 7 and + ! tr_prep_ns never matches a He5 species, so the THe35 product would + ! go to species slot 7 (past tr_pnf's (NSMAX,NNFMAX,NRMAX) arrays when + ! NSMAX < 7). Refused until a He5 species can be defined. + WRITE(6,*) 'XX Error: libnf: model_pnf=14 is not available: ', & + 'no He5 species can be defined' + STOP CASE DEFAULT - WRITE(6,*) 'XX Error libnb: undefined model_pnf: model_pnf=',model_pnf + WRITE(6,*) 'XX Error libnf: undefined model_pnf: model_pnf=',model_pnf STOP END SELECT @@ -236,15 +247,15 @@ SUBROUTINE set_usigmav_nf nsp_idnf(id_nf_dd2)=NS_H END IF eng_idnf(id_nf_dd2)=3.02D3*RKEV - enn_idnf(id_nf_dd3)=0.D0 + enn_idnf(id_nf_dd2)=0.D0 ns1_idnf(id_nf_dd3)=NS_D ns2_idnf(id_nf_dd3)=NS_D wgt_idnf(id_nf_dd3)=0.5D0 - IF(model_pnf.EQ.12) THEN - nsp_idnf(id_nf_dd3)=NS_He3 - ELSE + IF(model_pnf.EQ.12.OR.model_pnf.EQ.14) THEN nsp_idnf(id_nf_dd3)=NS_He4 + ELSE + nsp_idnf(id_nf_dd3)=NS_He3 END IF eng_idnf(id_nf_dd3)=0.82D3*RKEV enn_idnf(id_nf_dd3)=2.45D3*RKEV @@ -377,7 +388,7 @@ FUNCTION sigma_nf(id_nf,energy) RETURN END FUNCTION sigma_nf - ! --- reaction rate of nuclear fusion: sigmav --- + ! --- reaction rate of nuclear fusion: sigmav [m^3/s] --- ! --- as a function of temperature in keV FUNCTION sigmav_nf(id_nf,temperature) @@ -425,6 +436,12 @@ FUNCTION sigmav_nf(id_nf,temperature) WRITE(6,'(A,ES12.4)') ' temperature=',temperature STOP END IF + ! The svnf_* tables are in cm^3/s (D-T peaks at 8.7e-16 cm^3/s, as in + ! Bosch-Hale); tr_pnf needs m^3/s (its SNF has RN in 1e20 m^-3), and + ! the menu's B command prints the value. One conversion at the + ! computed exit covers all five tables; the low-temperature branch + ! returns 0 before it. + sigmav_nf=sigmav_nf*1.D-6 END FUNCTION sigmav_nf ! --- sigmav_nf by integeral over energy --- diff --git a/trx/trinit.f90 b/trx/trinit.f90 index b2ecde2de..64e4e2e2c 100644 --- a/trx/trinit.f90 +++ b/trx/trinit.f90 @@ -495,6 +495,7 @@ SUBROUTINE tr_init ! T + He3 -> He4 + p + n THe31 THe32 ! T + He3 -> He4 + D THe33 THe34 ! T + He3 -> He5 + p ! He4 + p THe35 THe36 + ! not available: no He5 species can be defined model_pnf=0 diff --git a/trx/trmenu.f90 b/trx/trmenu.f90 index 84e5df54d..9d3e6194c 100644 --- a/trx/trmenu.f90 +++ b/trx/trmenu.f90 @@ -154,7 +154,7 @@ SUBROUTINE tr_menu 7777 CONTINUE WRITE(6,*) '## INPUT id,T(kev)T' READ(5,*,END=7779) id,temp - WRITE(6,*) ' sigma=',sigmav_nf(id,temp) + WRITE(6,*) ' sigmav [m^3/s]=',sigmav_nf(id,temp) GOTO 7777 7779 CONTINUE diff --git a/trx/trpnf.f90 b/trx/trpnf.f90 index 6c49468d8..185eb7a91 100644 --- a/trx/trpnf.f90 +++ b/trx/trpnf.f90 @@ -66,6 +66,8 @@ SUBROUTINE tr_pnf PNFCL_NSNNFNR(1:NSMAX,1:NNFMAX,1:NRMAX)=0.D0 ! collisional transfer in SNFNN_NNFNR(1:NNFMAX,1:NRMAX)=0.D0 ! neutron number PNFNN_NNFNR(1:NNFMAX,1:NRMAX)=0.D0 ! neutron power + PNF_NSNNFNR(1:NSMAX,1:NNFMAX,1:NRMAX)=0.D0 ! fusion power + PNFIN_NSNNFNR(1:NSMAX,1:NNFMAX,1:NRMAX)=0.D0 ! fusion power in DO nnf=1,nnfmax id_nf=id_nf_nnf(nnf) @@ -105,10 +107,10 @@ SUBROUTINE tr_pnf DO nnf=1,nnfmax nsp=nsp_nnf(nnf) WF = RW(NR,NNBMAX+NNF) - VF =SQRT(2.D0*eng_nnf(nnf)*RKEV/(PA(ns)*AMP)) + VF =SQRT(2.D0*eng_nnf(nnf)/(PA(nsp)*AMP)) ! eng_nnf is in J already HYF=HY(VF/VCR) - TAUS = 0.2D0*PA(ns)*ABS(TE)**1.5D0 & - /(PZ(ns)**2*ANE*COULOG(1,ns,ANE,TE)) + TAUS = 0.2D0*PA(nsp)*ABS(TE)**1.5D0 & + /(PZ(nsp)**2*ANE*COULOG(1,nsp,ANE,TE)) TAUF(NNF,NR)= 0.5D0*TAUS*(1.D0-HYF) END DO END DO