summaryrefslogtreecommitdiff
path: root/Dragon/src/FLULPN.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/FLULPN.f')
-rw-r--r--Dragon/src/FLULPN.f227
1 files changed, 227 insertions, 0 deletions
diff --git a/Dragon/src/FLULPN.f b/Dragon/src/FLULPN.f
new file mode 100644
index 0000000..7e2275c
--- /dev/null
+++ b/Dragon/src/FLULPN.f
@@ -0,0 +1,227 @@
+*DECK FLULPN
+ SUBROUTINE FLULPN(IPMACR,NUNKNO,OPTION,TYPE,NGRP,NREG,NMAT,
+ 1 VOL,MATCOD,NMERG,IMERG,KEYFLX,FLUX,B2,IMPX,DIFHET,DHOM)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Calculation of heterogeneous leakage coefficients using the Todorova
+* approximation.
+*
+*Copyright:
+* Copyright (C) 2022 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
+* IPMACR pointer to the macrolib LCM object (L_MACROLIB signature).
+* NUNKNO number of flux/current unknowns.
+* OPTION type of leakage coefficients; can be 'P0' (P-0), 'P1' (P-1),
+* 'P0TR' (P-0 with transport correction).
+* TYPE type of buckling iteration.
+* Can be 'DIFF' (do a P0 calculation of DIFHET and exit);
+* other (do another type of calculation).
+* NGRP number of groups.
+* NREG number of volumes.
+* NMAT number of material mixtures.
+* VOL volumes.
+* MATCOD mixture number of each volume.
+* NMERG number of leakage zones.
+* IMERG leakage zone index in each material mixture zone.
+* KEYFLX position of each flux in the unknown vector.
+* FLUX direct unknown vector.
+* B2 buckling.
+* IMPX print flag.
+*
+*Parameters: output
+* DIFHET heterogeneous diffusion coefficients.
+* DHOM homogeneous diffusion coefficients.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ CHARACTER*4 OPTION,TYPE
+ TYPE(C_PTR) IPMACR
+ INTEGER NUNKNO,NGRP,NREG,NMAT,MATCOD(NREG),NMERG,IMERG(NMAT),
+ 1 KEYFLX(NREG),IMPX
+ REAL VOL(NREG),FLUX(NUNKNO,NGRP),B2,DIFHET(NMERG,NGRP),DHOM(NGRP)
+*----
+* LOCAL VARIABLES
+*----
+ TYPE(C_PTR) JPMACR,KPMACR
+ CHARACTER HSMG*131
+ DOUBLE PRECISION B1GAMA,DDELN1,DDELN2,DDELD1,B2HOM,ST2,STR,GAMMA
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS
+ REAL, ALLOCATABLE, DIMENSION(:) :: WORK
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: ST,FLXIN
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: SCAT1
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: STOD
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(FLXIN(NMAT,NGRP),ST(NMAT,NGRP),SCAT1(NMAT,NGRP,NGRP))
+*----
+* INITIALIZATION
+*----
+ IF((IMPX.GT.0).AND.(TYPE.EQ.'DIFF')) THEN
+ WRITE (6,100)
+ ELSE IF(IMPX.GT.0) THEN
+ WRITE (6,110) OPTION,TYPE,B2
+ ENDIF
+ ST(:NMAT,:NGRP)=0.0
+ SCAT1(:NMAT,:NGRP,:NGRP)=0.0
+*----
+* RECOVER MACROSCOPIC CROSS SECTIONS
+*----
+ ALLOCATE(IJJ(NMAT),NJJ(NMAT),IPOS(NMAT),WORK(NMAT*NGRP))
+ JPMACR=LCMGID(IPMACR,'GROUP')
+ DO IGR=1,NGRP
+ KPMACR=LCMGIL(JPMACR,IGR)
+ CALL LCMGET(KPMACR,'NTOT0',ST(1,IGR))
+ CALL LCMLEN(KPMACR,'SCAT01',ILONG,ITYLCM)
+ IF((ILONG.NE.0).AND.(OPTION.NE.'P0').AND.(OPTION.NE.'B0')) THEN
+ CALL LCMGET(KPMACR,'IJJS01',IJJ)
+ CALL LCMGET(KPMACR,'NJJS01',NJJ)
+ CALL LCMGET(KPMACR,'IPOS01',IPOS)
+ CALL LCMGET(KPMACR,'SCAT01',WORK)
+ DO IBM=1,NMAT
+ IPO=IPOS(IBM)
+ DO JGR=IJJ(IBM),IJJ(IBM)-NJJ(IBM)+1,-1
+ SCAT1(IBM,IGR,JGR)=WORK(IPO) ! IGR <-- JGR
+ IPO=IPO+1
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ DEALLOCATE(WORK,IPOS,NJJ,IJJ)
+*----
+* RECOVER INTEGRATED FLUX
+*----
+ FLXIN(:NMAT,:NGRP)=0.0
+ DO IGR=1,NGRP
+ DO IBM=1,NMAT
+ DO I=1,NREG
+ IND=KEYFLX(I)
+ IF((MATCOD(I).EQ.IBM).AND.(IND.GT.0)) THEN
+ FLXIN(IBM,IGR)=FLXIN(IBM,IGR)+FLUX(IND,IGR)*VOL(I)
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ IF((OPTION.EQ.'LKRD').OR.(OPTION.EQ.'RHS')) GO TO 10
+*----
+* MAIN LOOP OVER LEAKAGE ZONES
+*----
+ B2HOM=DBLE(B2)
+ GAMMA=1.0D0
+ DO INM=1,NMERG
+ IF((OPTION.EQ.'P0').OR.(OPTION.EQ.'B0')) THEN
+* P0 or B0 approximation
+ DO IGR=1,NGRP
+ DDELN1=0.D0
+ DDELD1=0.D0
+ DO IBM=1,NMAT
+ IF(IMERG(IBM).EQ.INM) THEN
+ DDELN1=DDELN1+ST(IBM,IGR)*FLXIN(IBM,IGR)
+ DDELD1=DDELD1+FLXIN(IBM,IGR)
+ ENDIF
+ ENDDO
+ ST2=DDELN1/DDELD1
+ IF(OPTION.EQ.'B0') GAMMA=B1GAMA(2,B2HOM,ST2)
+ DIFHET(INM,IGR)=REAL(1.0D0/(3.0D0*GAMMA*ST2))
+ ENDDO
+ ELSE IF((OPTION.EQ.'P0TR').OR.(OPTION.EQ.'B0TR').OR.
+ 1 (TYPE.EQ.'DIFF')) THEN
+* Outscatter approximation
+ DO IGR=1,NGRP
+ DDELN1=0.D0
+ DDELN2=0.D0
+ DDELD1=0.D0
+ DO IBM=1,NMAT
+ IF(IMERG(IBM).EQ.INM) THEN
+ DDELN1=DDELN1+ST(IBM,IGR)*FLXIN(IBM,IGR)
+ DO JGR=1,NGRP
+ DDELN2=DDELN2+SCAT1(IBM,JGR,IGR)*FLXIN(IBM,IGR)
+ ENDDO
+ DDELD1=DDELD1+FLXIN(IBM,IGR)
+ ENDIF
+ ENDDO
+ ST2=DDELN1/DDELD1
+ IF(OPTION.EQ.'B0TR') GAMMA=B1GAMA(2,B2HOM,ST2)
+ STR=(GAMMA*DDELN1-DDELN2)/DDELD1
+ DIFHET(INM,IGR)=REAL(1.0D0/(3.0D0*STR))
+ ENDDO
+ ELSE IF((OPTION.EQ.'P1').OR.(OPTION.EQ.'B1')) THEN
+* Inscatter approximation
+ ALLOCATE(STOD(NGRP,NGRP+1))
+ STOD(:NGRP,:NGRP+1)=0.0D0
+ DO IGR=1,NGRP
+ IF(OPTION.EQ.'B1') THEN
+ DDELN1=0.D0
+ DDELD1=0.D0
+ DO IBM=1,NMAT
+ IF(IMERG(IBM).EQ.INM) THEN
+ DDELN1=DDELN1+ST(IBM,IGR)*FLXIN(IBM,IGR)
+ DDELD1=DDELD1+FLXIN(IBM,IGR)
+ ENDIF
+ ENDDO
+ ST2=DDELN1/DDELD1
+ GAMMA=B1GAMA(2,B2HOM,ST2)
+ ENDIF
+ DO IBM=1,NMAT
+ IF(IMERG(IBM).EQ.INM) THEN
+ STOD(IGR,IGR)=STOD(IGR,IGR)+GAMMA*ST(IBM,IGR)*
+ 1 FLXIN(IBM,IGR)
+ DO JGR=1,NGRP
+ STOD(IGR,JGR)=STOD(IGR,JGR)-SCAT1(IBM,IGR,JGR)*
+ 1 FLXIN(IBM,JGR)
+ ENDDO
+ STOD(IGR,NGRP+1)=STOD(IGR,NGRP+1)+FLXIN(IBM,IGR)/3.0D0
+ ENDIF
+ ENDDO
+ ENDDO
+ CALL ALSBD(NGRP,1,STOD,IER,NGRP)
+ IF(IER.NE.0) CALL XABORT('FLULPN: SINGULAR MATRIX.')
+ DO IGR=1,NGRP
+ DIFHET(INM,IGR)=REAL(STOD(IGR,NGRP+1))
+ ENDDO
+ DEALLOCATE(STOD)
+ ELSE
+ WRITE(HSMG,'(15HFLULPN: OPTION ,A,23H IS INVALID WITH TODORO,
+ 1 17HVA APPROXIMATION.)') OPTION
+ CALL XABORT(HSMG)
+ ENDIF
+ ENDDO
+*----
+* COMPUTE THE HOMOGENEOUS LEAKAGE COEFFICIENTS
+*----
+ 10 DO IGR=1,NGRP
+ DHOM(IGR)=0.0
+ FLTOT=0.0
+ DO IBM=1,NMAT
+ INM=IMERG(IBM)
+ DHOM(IGR)=DHOM(IGR)+FLXIN(IBM,IGR)*DIFHET(INM,IGR)
+ FLTOT=FLTOT+FLXIN(IBM,IGR)
+ ENDDO
+ DHOM(IGR)=DHOM(IGR)/FLTOT
+ ENDDO
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(SCAT1,ST,FLXIN)
+ RETURN
+*
+ 100 FORMAT(/54H FLULPN: OUTSCATTER DIFFUSION COEFFICIENT CALCULATION.)
+ 110 FORMAT(/21H FLULPN: SOLUTION OF ,A4,21H EQUATIONS WITH TYPE ,A4,
+ 1 10H BUCKLING=,1P,E12.4)
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 15:34:29 +0000