Merge pull request #275 from PrincetonUniversity/feature/b3d_fusion_p…

…ower BEAMS3D: Added fusion power integrals to screen output.
PrincetonUniversity · Aug 9, 2024 · 0c123c3 · 0c123c3
2 parents 7ce89b5 + 84f41e3
commit 0c123c3
Showing 1 changed file with 22 additions and 5 deletions.
diff --git a/BEAMS3D/Sources/beams3d_init_fusion.f90 b/BEAMS3D/Sources/beams3d_init_fusion.f90
@@ -35,7 +35,7 @@ SUBROUTINE beams3d_init_fusion
 !-----------------------------------------------------------------------
       IMPLICIT NONE
       LOGICAL :: lfusion_DHe3
-      INTEGER :: s,i,j,k,k1,k2, nr1, nphi1, nz1, l
+      INTEGER :: s,i,j,k,k1,k2, nr1, nphi1, nz1, l, nfp
       INTEGER :: minik(2)
       INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: n3d
       REAL(rprec) :: maxrateDT, maxrateDDT, maxrateDDHe, maxrateDHe3, &
@@ -92,6 +92,7 @@ SUBROUTINE beams3d_init_fusion
       nr1   = nr - 1
       nphi1 = nphi - 1
       nz1   = nz - 1
+      nfp = NINT(pi2/phiaxis(nphi))
 
       ! We need to first define the reaction rate over the grid
       CALL mpialloc(l3d,      nr1, nphi1, nz1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_l3d)
@@ -196,10 +197,26 @@ SUBROUTINE beams3d_init_fusion
       ! Output reaction rate info
       IF (myworkid == master) THEN
          IF (lverb) THEN
-            IF (lfusion_alpha) WRITE(6, '(A,ES11.4,A)') '      max D + T -> He4   rate: ', maxrateDT,' [part/(m^3 s)]'
-            IF (lfusion_tritium .or. lfusion_proton) WRITE(6, '(A,ES11.4,A)') '      max D + D -> T + p rate: ', maxrateDDT,' [part/(m^3 s)]'
-            IF (lfusion_He3) WRITE(6, '(A,ES11.4,A)')  '      max D + D -> He3   rate: ', maxrateDDHe,' [part/(m^3 s)]'
-            IF (lfusion_DHe3) WRITE(6, '(A,ES11.4,A)') '      max D +He3-> He4 + p   rate: ', maxrateDHe3,' [part/(m^3 s)]'
+            IF (lfusion_alpha) THEN
+               WRITE(6, '(A,ES11.4,A)') '      max D + T -> He4   rate: ', maxrateDT,' [part/(m^3 s)]'
+               WRITE(6, '(A,F7.2,A)')   '          He4(3.52 MeV) Power: ', SUM(rateDT)*(3.52)*e_charge*nfp,' [MW]'
+               WRITE(6, '(A,F7.2,A)')   '           n(14.06 MeV) Power: ', SUM(rateDT)*(14.06)*e_charge*nfp,' [MW]'
+            END IF
+            IF (lfusion_tritium .or. lfusion_proton) THEN
+               WRITE(6, '(A,ES11.4,A)') '      max D + D -> T + p rate: ', maxrateDDT,' [part/(m^3 s)]'
+               WRITE(6, '(A,F7.2,A)')   '            T(1.01 MeV) Power: ', SUM(rateDDT)*(1.01)*e_charge*nfp,' [MW]'
+               WRITE(6, '(A,F7.2,A)')   '            p(3.02 MeV) Power: ', SUM(rateDDT)*(3.02)*e_charge*nfp,' [MW]'
+            END IF
+            IF (lfusion_He3) THEN
+               WRITE(6, '(A,ES11.4,A)') '      max D + D -> He3   rate: ', maxrateDDHe,' [part/(m^3 s)]'
+               WRITE(6, '(A,F7.2,A)')   '          He3(0.82 MeV) Power: ', SUM(rateDDHe)*(0.82)*e_charge*nfp,' [MW]'
+               WRITE(6, '(A,F7.2,A)')   '            n(2.45 MeV) Power: ', SUM(rateDDHe)*(2.45)*e_charge*nfp,' [MW]'
+            END IF
+            IF (lfusion_DHe3) THEN
+               WRITE(6, '(A,ES11.4,A)') '      max D + He3 -> He4 + p   rate: ', maxrateDHe3,' [part/(m^3 s)]'
+               WRITE(6, '(A,F7.2,A)')   '                 He4(3.6 MeV) Power: ', SUM(rateDHe3)*(3.60)*e_charge*nfp,' [MW]'
+               WRITE(6, '(A,F7.2,A)')   '                  p(14.7 MeV) Power: ', SUM(rateDHe3)*(14.7)*e_charge*nfp,' [MW]'
+            END IF
          END IF
       END IF