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/FLULPN.f | 227 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 Dragon/src/FLULPN.f (limited to 'Dragon/src/FLULPN.f') 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 -- cgit v1.2.3