summaryrefslogtreecommitdiff
path: root/Dragon/src/MUSA.f90
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/MUSA.f90')
-rw-r--r--Dragon/src/MUSA.f90287
1 files changed, 287 insertions, 0 deletions
diff --git a/Dragon/src/MUSA.f90 b/Dragon/src/MUSA.f90
new file mode 100644
index 0000000..aa650dd
--- /dev/null
+++ b/Dragon/src/MUSA.f90
@@ -0,0 +1,287 @@
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! Calculation of cellwise scattering-reduced collision, escape and
+! transmission probabilities for the current iteration method in the
+! multicell surfacic approximation.
+!
+!Copyright:
+! Copyright (C) 2025 Ecole Polytechnique de Montreal
+! This library is free software; you can redistribute it and/or
+! modify it under the terms of the GNU Lesser General Public
+! License as published by the Free Software Foundation; either
+! version 2.1 of the License, or (at your option) any later version
+!
+!Author(s): A. Hebert
+!
+!Parameters: input
+! IPSYS pointer to the system matrices.
+! IPTRK pointer to the tracking (L_TRACK signature).
+! IFTRAK tracking file unit.
+! IMPX print flag (equal to zero for no print).
+! NREG total number of merged regions for which specific values
+! of the neutron flux and reactions rates are required.
+! NBMIX number of mixtures.
+! SIGT0 total macroscopic cross sections ordered by mixture.
+! SIGW0 within-group scattering macroscopic cross section ordered
+! by mixture.
+! NBATCH number of tracks dispached in eack OpenMP core.
+! TITREC title.
+! NALBP number of multigroup physical albedos.
+! ALBP multigroup physical albedos.
+!
+!-----------------------------------------------------------------------
+!
+SUBROUTINE MUSA(IPSYS,IPTRK,IFTRAK,IMPX,NREG,NBMIX,SIGT0,SIGW0,NBATCH, &
+ & TITREC,NALBP,ALBP)
+ USE GANLIB
+ !----
+ ! SUBROUTINE ARGUMENTS
+ !----
+ TYPE(C_PTR) IPSYS,IPTRK
+ INTEGER IFTRAK,IMPX,NREG,NBMIX,NBATCH,NALBP
+ REAL SIGT0(0:NBMIX),SIGW0(0:NBMIX),ALBP(NALBP)
+ CHARACTER TITREC*72
+ !----
+ ! LOCAL VARIABLES
+ !----
+ PARAMETER (EPS1=1.0E-4,NSTATE=40)
+ TYPE(C_PTR) JPTRK,KPTRK
+ INTEGER ISTATT(NSTATE),NNPSYS(1)
+ CHARACTER TITRE2*72
+ logical LSKIP
+ !----
+ ! ALLOCATABLE ARRAYS
+ !----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: MATALB,NMC_NODE,NMC_SURF,MAT2,IGEN,INUM
+ REAL, ALLOCATABLE, DIMENSION(:) :: SIGT2,SIGW2,WORK
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: VOLSUR
+ REAL, POINTER, DIMENSION(:) :: PSSW,PSJW,PISW,PIJW
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DPROB,DPROBX
+ TYPE(C_PTR) :: PSSW_PTR,PSJW_PTR,PISW_PTR,PIJW_PTR
+ !
+ IND(I,J) = MAX(I+J3+1,J+J3+1)*(MAX(I+J3+1,J+J3+1)-1)/2 &
+ & + MIN(I+J3+1,J+J3+1)
+ !
+ WPR(I,J)= REAL(DPROB( IND(I,J),1 ) / DPROB( IND(I,0),1 ))
+ !----
+ ! BICKLEY FLAG
+ !----
+ SAVE IBICKL
+ DATA IBICKL/0/
+ !----
+ ! RECOVER BICKLEY TABLES
+ !----
+ IF(IBICKL.EQ.0) THEN
+ CALL XDRTA2
+ IBICKL=1
+ ENDIF
+ !----
+ ! RECOVER SALT SPECIFIC PARAMETERS
+ !----
+ REWIND IFTRAK
+ CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATT)
+ IF(NREG.NE.ISTATT(1)) THEN
+ CALL XABORT('MUSA: STATE VECTOR HAS INVALID # OF ZONES.')
+ ENDIF
+ NMACRO=ISTATT(24) ! NGEN
+ NMCEL=NMACRO
+ NMERGE=NMACRO
+ NGEN=NMACRO
+ ALLOCATE(IGEN(NMERGE),INUM(NMCEL))
+ DO IK=1,NMERGE
+ IGEN(IK)=IK
+ ENDDO
+ DO IK=1,NMCEL
+ INUM(IK)=IK
+ ENDDO
+ IF(NMACRO.EQ.0) CALL XABORT('MUSA: MUST MODULE TRACKING IS MANDATORY.')
+ ALLOCATE(NMC_NODE(NMACRO+1),NMC_SURF(NMACRO+1))
+ JPTRK=LCMGID(IPTRK,'MACRO-TRACK')
+ CALL LCMGET(IPTRK,'NMC_NODE',NMC_NODE)
+ CALL LCMGET(IPTRK,'NMC_SURF',NMC_SURF)
+ NMIX=NMC_SURF(NMACRO+1)
+ NIFR=NMC_SURF(NMACRO+1)
+ !----
+ ! LOOP OVER MACRO GEOMETRIES AND COMPUTE PIJ MATRICES USING EXCELP
+ !----
+ J1=0
+ NMIX=0
+ NPIJ=0
+ NPIS=0
+ NPSS=0
+ DO IMACRO=1,NMACRO
+ J2=NMC_NODE(IMACRO+1)-NMC_NODE(IMACRO)
+ J3=NMC_SURF(IMACRO+1)-NMC_SURF(IMACRO)
+ J1=J1+J2
+ NMIX=NMIX+J3
+ NPIJ=NPIJ+J2*J2
+ NPIS=NPIS+J2*J3
+ NPSS=NPSS+J3*J3
+ ENDDO
+ IF(J1.NE.NREG) CALL XABORT('MUSA: INVALID NREG.')
+ IF(NMIX.NE.NMC_SURF(NMACRO+1)) CALL XABORT('MUSA: INVALID NMIX.')
+ PIJW_PTR=LCMARA(NPIJ)
+ PISW_PTR=LCMARA(NPIS)
+ PSJW_PTR=LCMARA(NPIS)
+ PSSW_PTR=LCMARA(NPSS)
+ CALL C_F_POINTER(PIJW_PTR,PIJW,(/ NPIJ /))
+ CALL C_F_POINTER(PISW_PTR,PISW,(/ NPIS /))
+ CALL C_F_POINTER(PSJW_PTR,PSJW,(/ NPIS /))
+ CALL C_F_POINTER(PSSW_PTR,PSSW,(/ NPSS /))
+ J1=0
+ IPIJ=0
+ IPIS=0
+ IPSS=0
+ DO IMACRO=1,NMACRO
+ J2=NMC_NODE(IMACRO+1)-NMC_NODE(IMACRO)
+ J3=NMC_SURF(IMACRO+1)-NMC_SURF(IMACRO)
+ N2PRO=(J2+J3+1)**2
+ WRITE(TITRE2,'(A,9H -- MACRO,I5.5)') TRIM(TITREC),IMACRO
+ KPTRK=LCMGIL(JPTRK,IMACRO)
+ KNORM=4 ! use HELIOS-type normalization
+ NNPSYS(1)=1
+ ALLOCATE(MAT2(J2),SIGT2(J2),SIGW2(J2))
+ ALLOCATE(MATALB(-J3:J2),VOLSUR(-J3:J2,1),DPROB(N2PRO,1),DPROBX(N2PRO,1))
+ CALL LCMGET(KPTRK,'MATCOD',MAT2)
+ CALL EXCELP(KPTRK,IFTRAK,IMPX,J3,J2,NBMIX,MAT2,KNORM,SIGT0,1,N2PRO, &
+ & 1,NNPSYS(1),NBATCH,TITRE2,NALBP,ALBP,MATALB,VOLSUR,DPROB,DPROBX)
+ !----
+ ! CHECK IF SCATTERING REDUCTION IS REQUIRED
+ !----
+ DO I=1,J2
+ SIGT2(I)=SIGT0(MAT2(I)) ! sigt by node
+ SIGW2(I)=SIGW0(MAT2(I)) ! sigw by node
+ ENDDO
+ LSKIP=.TRUE.
+ DO I=1,J2
+ LSKIP=LSKIP.AND.(SIGW2(I).EQ.0.0)
+ ENDDO
+ !----
+ ! SCATTERING REDUCTION IF LSKIP=.FALSE.
+ !----
+ IF(LSKIP) THEN
+ ! DO NOT PERFORM SCATTERING REDUCTION.
+ DO I=1,J2
+ DO J=1,J2
+ IF(SIGT2(J).EQ.0.0) THEN
+ PIJW(IPIJ+(J-1)*J2+I)=WPR(I,J)
+ ELSE
+ PIJW(IPIJ+(J-1)*J2+I)=WPR(I,J)/SIGT2(J)
+ ENDIF
+ ENDDO
+ ENDDO
+ DO I=1,J2
+ DO JC=1,J3
+ PISW(IPIS+(JC-1)*J2+I)=WPR(I,-JC)
+ IF(SIGT2(I).EQ.0.0) THEN
+ PSJW(IPIS+(I-1)*J3+JC)=WPR(-JC,I)
+ ELSE
+ PSJW(IPIS+(I-1)*J3+JC)=WPR(-JC,I)/SIGT2(I)
+ ENDIF
+ ENDDO
+ ENDDO
+ DO IC=1,J3
+ DO JC=1,J3
+ PSSW(IPSS+(JC-1)*J3+IC)=WPR(-IC,-JC)
+ ENDDO
+ ENDDO
+ ELSE
+ ! COMPUTE THE SCATTERING-REDUCED COLLISION AND ESCAPE MATRICES.
+ DO I=1,J2
+ DO J=1,J2
+ IF(SIGT2(J).EQ.0.0) THEN
+ PIJW(IPIJ+(J-1)*J2+I)=0.0
+ ELSE
+ PIJW(IPIJ+(J-1)*J2+I)=-WPR(I,J)*SIGW2(J)/SIGT2(J)
+ ENDIF
+ ENDDO
+ PIJW(IPIJ+(I-1)*J2+I)=1.0+PIJW(IPIJ+(I-1)*J2+I)
+ ENDDO
+ CALL ALINV(J2,PIJW(IPIJ+1),J2,IER)
+ IF(IER.NE.0) CALL XABORT('MUSA: SINGULAR MATRIX.')
+ ALLOCATE(WORK(J2))
+ DO I=1,J2
+ DO K=1,J2
+ WORK(K)=PIJW(IPIJ+(K-1)*J2+I)
+ ENDDO
+ DO J=1,J2
+ WGAR=0.0
+ DO K=1,J2
+ IF(SIGT2(J).EQ.0.0) THEN
+ WGAR=WGAR+WORK(K)*WPR(K,J)
+ ELSE
+ WGAR=WGAR+WORK(K)*WPR(K,J)/SIGT2(J)
+ ENDIF
+ ENDDO
+ PIJW(IPIJ+(J-1)*J2+I)=WGAR
+ ENDDO
+ DO JC=1,J3
+ WGAR=0.0
+ DO K=1,J2
+ WGAR=WGAR+WORK(K)*WPR(K,-JC)
+ ENDDO
+ PISW(IPIS+(JC-1)*J2+I)=WGAR
+ ENDDO
+ ENDDO
+ DEALLOCATE(WORK)
+ !
+ ! COMPUTE THE SCATTERING-REDUCED COLLISION PROBABILITY MATRIX
+ ! FOR INCOMING NEUTRONS.
+ DO IC=1,J3
+ DO J=1,J2
+ IF(SIGT2(J).EQ.0.0) THEN
+ WGAR=WPR(-IC,J)
+ ELSE
+ WGAR=WPR(-IC,J)/SIGT2(J)
+ ENDIF
+ DO K=1,J2
+ IF(SIGT2(K).NE.0.0) THEN
+ WGAR=WGAR+WPR(-IC,K)*PIJW(IPIJ+(J-1)*J2+K)*SIGW2(K)/SIGT2(K)
+ ENDIF
+ ENDDO
+ PSJW(IPIS+(J-1)*J3+IC)=WGAR
+ ENDDO
+ ENDDO
+ !
+ ! COMPUTE THE SCATTERING-REDUCED TRANSMISSION PROBABILITY MATRIX.
+ DO IC=1,J3
+ DO JC=1,J3
+ WGAR=WPR(-IC,-JC)
+ DO K=1,J2
+ IF(SIGT2(K).NE.0.0) THEN
+ WGAR=WGAR+WPR(-IC,K)*PISW(IPIS+(JC-1)*J2+K)*SIGW2(K)/SIGT2(K)
+ ENDIF
+ ENDDO
+ PSSW(IPSS+(JC-1)*J3+IC)=WGAR
+ ENDDO
+ ENDDO
+ ENDIF
+ DEALLOCATE(DPROBX,DPROB,VOLSUR,MATALB)
+ IF(IMPX.GE.8) THEN
+ IF(LSKIP) THEN
+ IN=1
+ ELSE
+ IN=2
+ ENDIF
+ CALL SYBPRX(IN,J3,J2,IMACRO,SIGT2,SIGW2,PIJW(IPIJ+1),PISW(IPIS+1), &
+ & PSJW(IPIS+1),PSSW(IPSS+1))
+ ENDIF
+ DEALLOCATE(SIGW2,SIGT2,MAT2)
+ J1=J1+J2
+ IPIJ=IPIJ+J2*J2
+ IPIS=IPIS+J2*J3
+ IPSS=IPSS+J3*J3
+ ENDDO
+ ! end of SYB004 equivalent
+ CALL LCMPPD(IPSYS,'PSSW$SYBIL',NPSS,2,PSSW_PTR)
+ CALL LCMPPD(IPSYS,'PSJW$SYBIL',NPIS,2,PSJW_PTR)
+ CALL LCMPPD(IPSYS,'PISW$SYBIL',NPIS,2,PISW_PTR)
+ CALL LCMPPD(IPSYS,'PIJW$SYBIL',NPIJ,2,PIJW_PTR)
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(/31H MUSA: PIJ INFORMATION IN GROUP)')
+ CALL LCMLIB(IPSYS)
+ ENDIF
+ RETURN
+END SUBROUTINE MUSA
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 18:07:01 +0000