From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/MUSA.f90 | 287 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 287 insertions(+) create mode 100644 Dragon/src/MUSA.f90 (limited to 'Dragon/src/MUSA.f90') 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 -- cgit v1.2.3