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/PIJXL3.f | 449 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 449 insertions(+) create mode 100644 Dragon/src/PIJXL3.f (limited to 'Dragon/src/PIJXL3.f') diff --git a/Dragon/src/PIJXL3.f b/Dragon/src/PIJXL3.f new file mode 100644 index 0000000..1d79b45 --- /dev/null +++ b/Dragon/src/PIJXL3.f @@ -0,0 +1,449 @@ +*DECK PIJXL3 + SUBROUTINE PIJXL3(IPTRK,IPRT,NGRP,NANI,NBMIX,NPSYS,NRENOR,LEAKSW, + > XSSIGT,XSSIGW,NELPIJ,PIJ) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Calculation of the collision probabilities in EXCELL without producing +* a tracking file. Based on subroutine XL3TRK in DRAGON 3.4. +* +*Copyright: +* Copyright (C) 2005 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): R. Roy and A. Hebert +* +*Parameters: input +* IPTRK pointer to the tracking (L_TRACK signature). +* IPRT print flag (equal to zero for no print). +* NGRP number of energy groups. +* NANI number of Legendre orders (usually equal to one). +* NBMIX number of mixtures. +* NPSYS index set to zero if a group is not to be processed. Usually, +* NPSYS(I)=I. +* NRENOR normalization scheme for PIJ matrices. +* LEAKSW leakage flag (=.true. if neutron leakage through external +* boundary is present). +* XSSIGT total macroscopic cross sections ordered by mixture. +* XSSIGW P0 within-group scattering macroscopic cross sections +* ordered by mixture. +* NELPIJ number of elements in symmetrized pij matrix. +* +*Parameters: output +* PIJ reduced and symmetrized collision probabilities. +* +*----------------------------------------------------------------------- +* + USE GANLIB + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + LOGICAL LEAKSW + TYPE(C_PTR) IPTRK + INTEGER IPRT,NGRP,NANI,NBMIX,NPSYS(NGRP),NRENOR,NELPIJ + REAL XSSIGT(0:NBMIX,NGRP),XSSIGW(0:NBMIX,NANI,NGRP), + > PIJ(NELPIJ,NGRP) +*---- +* LOCAL VARIABLES +*---- + INTEGER IOUT,NALB,NSTATE,ICPALL,ICPEND + PARAMETER (IOUT=6,NALB=6,NSTATE=40,ICPALL=4,ICPEND=3) + INTEGER NDIM ,ISPEC ,NANGLE,NANGL ,ISYMM, NORE + INTEGER NALBG ,NC ,NTR ,NTZ ,ICL ,NSOUT , + > ITGEO ,NRMV ,NTY ,LTRK ,LINMAX,NUNK , + > MAXR ,NEXTGE,NTX ,NCOR ,NSUR ,NTOTCL, + > NVOL ,NV ,NS ,IGRP ,ISOUT ,ILONG , + > ITYPE ,INDPIJ,IIN ,IBM ,I ,J , + > NPIJ ,NREG ,NUNKMR + INTEGER ISTATE(NSTATE),LCLSYM(3) + INTEGER MXANGL,ICODE(NALB) + LOGICAL SWVOID,SWNZBC,LSKIP + REAL ALBOLD(NALB),EXTKOP(NSTATE),DENUSR, RCUTOF, + > CUTOFX, FACT +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: MINDIM,MAXDIM,ICORD,INDEL, + > KEYMRG,MATALB,MATMRG,ICUR,INCR,NUMERO,MATRT + REAL, ALLOCATABLE, DIMENSION(:) :: REMESH,VOLSUR,VOLMRG,CONV, + > TRKBEG,TRKDIR,FFACT,ANGLES,DENSTY + REAL, ALLOCATABLE, DIMENSION(:,:) :: SIGVOL,SIGTAL + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: LENGHT,VOLTRK, + > PSST,PSVT + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DBLPIJ,PCSCT +*---- +* INTRINSIC FUNCTION FOR POSITION IN CONDENSE PIJ MATRIX +*---- + INTEGER INDPOS + INDPOS(I,J)=MAX(I,J)*(MAX(I,J)-1)/2+MIN(I,J) +*---- +* READ THE GEOMETRY INFORMATION STORED ON IPTRK +*---- + ISTATE(:NSTATE)=0 + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + NORE=ISTATE(8) + LTRK=ISTATE(9)+1 + NANGLE=ISTATE(11) + ISYMM=ISTATE(12) + CALL LCMGET(IPTRK,'EXCELTRACKOP',EXTKOP) + CUTOFX=EXTKOP(1) + DENUSR=EXTKOP(2) + RCUTOF=EXTKOP(3) + CALL LCMSIX(IPTRK,'EXCELL',1) + ISTATE(:NSTATE)=0 + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + NDIM =ISTATE(1) + NSUR =-ISTATE(2) + NVOL =ISTATE(3) + NTOTCL =ISTATE(4) + MAXR =ISTATE(5) + NUNK =ISTATE(6) + NEXTGE =ISTATE(7) +*---- +* Intrinsic symmetries used in geometry +* Use these to simplify tracking unless +* NOSYMM tracking option activated +*---- + LCLSYM(1) =ISTATE(8) + LCLSYM(2) =ISTATE(9) + LCLSYM(3) =ISTATE(10) + IF(ISYMM .NE. 0) THEN + ISYMM=0 + IF(NDIM .EQ. 2) THEN + IF(LCLSYM(1) .NE. 0) THEN +*---- +* X SYMMETRY +*---- + ISYMM=2 + ENDIF + IF(LCLSYM(2) .NE. 0) THEN + IF(ISYMM .EQ. 0) THEN +*---- +* Y SYMMETRY +*---- + ISYMM=4 + ELSE +*---- +* X AND Y SYMMETRY +*---- + ISYMM=8 + ENDIF + ENDIF +*---- +* X-Y DIAGONAL SYMMETRY +*---- + ELSE + IF(LCLSYM(1) .NE. 0) THEN +*---- +* X SYMMETRY +*---- + ISYMM=2 + ENDIF + IF(LCLSYM(2) .NE. 0) THEN + IF(ISYMM .EQ. 0) THEN +*---- +* Y SYMMETRY +*---- + ISYMM=4 + ELSE +*---- +* X AND Y SYMMETRY +*---- + ISYMM=8 + ENDIF + ENDIF + IF(LCLSYM(3) .NE. 0) THEN +*---- +* Z SYMMETRY +*---- + ISYMM=ISYMM+16 + ENDIF + ENDIF + IF(ISYMM .EQ. 0) ISYMM=1 + ENDIF + ALLOCATE(MINDIM(NTOTCL),MAXDIM(NTOTCL),ICORD(NTOTCL), + > INDEL(4*NUNK),KEYMRG(NUNK),MATALB(NUNK)) + ALLOCATE(REMESH(MAXR),VOLSUR(NUNK)) + CALL LCMGET(IPTRK,'MINDIM ',MINDIM) + CALL LCMGET(IPTRK,'MAXDIM ',MAXDIM) + CALL LCMGET(IPTRK,'ICORD ',ICORD ) + CALL LCMGET(IPTRK,'INDEX ',INDEL ) + CALL LCMGET(IPTRK,'REMESH ',REMESH) + CALL LCMGET(IPTRK,'KEYMRG ',KEYMRG) + CALL LCMGET(IPTRK,'MATALB ',MATALB) + CALL LCMGET(IPTRK,'VOLSUR ',VOLSUR) + CALL LCMSIX(IPTRK,'EXCELL ',2) + CALL LCMGET(IPTRK,'ICODE ',ICODE ) + CALL LCMGET(IPTRK,'ALBEDO ',ALBOLD) +*---- +* VERIFY SYMMETRY AND STUDY TRACKING PARAMETERS. ARE THEY BASICALLY +* POSSIBLE ? +*---- + MXANGL=0 + IF(LTRK .EQ. 1)THEN + NCOR= 1 + IF(NDIM .EQ. 2) THEN + MXANGL=NANGLE + IF(ISYMM .GE. 2) THEN + NANGL = (NANGLE+1)/2 + ELSE + NANGL = NANGLE + ENDIF + IF( RCUTOF.GT.0.0 ) NCOR= 2 + ELSE IF(NDIM .EQ. 3) THEN + IF(MOD(NANGLE,2) .EQ. 1)THEN + NANGLE=NANGLE+1 + WRITE(IOUT,'(/31H MESS = ONLY EVEN # EQN ANGLES )') + ENDIF + IF(NANGLE .GT. 16)THEN + NANGLE=16 + WRITE(IOUT,'(/31H MESS = 16 IS MAX # EQN ANGLES )') + ENDIF + MXANGL=(NANGLE * (NANGLE+2)) / 2 + IF(NEXTGE .EQ. 1) THEN + NANGL = (NANGLE * (NANGLE+2)) / 8 + ELSE + IF(ISYMM .EQ. 8 .OR. ISYMM .EQ. 24) THEN + NANGL = (NANGLE * (NANGLE+2)) / 8 + ELSE IF(ISYMM .EQ. 2 .OR. ISYMM .EQ. 4 .OR. + > ISYMM .EQ. 18 .OR. ISYMM .EQ. 20 ) THEN + NANGL = (NANGLE * (NANGLE+2)) / 4 + ELSE + NANGL = (NANGLE * (NANGLE+2)) / 2 + ENDIF + ENDIF + IF(RCUTOF .GT. 0.0) NCOR= 4 + ENDIF + ELSEIF( LTRK.EQ.2 )THEN + NCOR = 1 + MXANGL=NANGLE + IF( NDIM.EQ.2 )THEN + NANGL = NANGLE + ELSEIF( NDIM.EQ.3 )THEN + CALL XABORT('PIJXL3: *TSPC* NOT AVAILABLE FOR 3-D GEOMETRY') + ENDIF + CUTOFX= RCUTOF + ENDIF + IF( IPRT.GE.1 ) THEN + WRITE(IOUT,6002) NANGL,ISYMM,CUTOFX,DENUSR,RCUTOF + ENDIF + IF( IPRT.GT.1 .AND. NEXTGE.EQ.0 )THEN +* +* IF PRINT REQUIRED AND OVERALL CARTESIAN GEOMETRY +* PRINT CARTESIAN REGION MAP + NTX= MAXDIM(1)-MINDIM(1) + NTY= MAXDIM(2)-MINDIM(2) + NTZ= MAXDIM(3)-MINDIM(3) + NTR=0 + DO 103 ICL=4,NTOTCL + NTR= MAX(NTR,MAXDIM(ICL)-MINDIM(ICL)+1) + 103 CONTINUE + CALL XELGPR(NDIM,NTX,NTY,NTZ,NTR,ISYMM, + > NSUR,NVOL,NTOTCL,MINDIM,MAXDIM, + > KEYMRG,INDEL,MATALB) + ENDIF + ALLOCATE(VOLTRK((NANGL+1)*NUNK)) +* + NV= NVOL + NS= -NSUR + ALLOCATE(VOLMRG(NUNK),MATMRG(NUNK)) + ITGEO=3 + CALL XELCMP( NS, NV, VOLSUR, MATALB, KEYMRG, + > NSOUT, NREG, VOLMRG, MATMRG, ITGEO,ICODE) + NUNKMR= NREG+NSOUT+1 + NPIJ= (NUNKMR*(NUNKMR+1))/2 + IF( IPRT .GT. 1 ) WRITE(IOUT,6000) (NGRP*NPIJ/128) + ALLOCATE(DBLPIJ(NPIJ,NGRP)) + IF( IPRT .GT. 1 ) WRITE(IOUT,6001) +* +* ALLOCATE AND CHARGE TOTAL XS PER REGION + ALLOCATE(SIGTAL(NUNKMR,NGRP),SIGVOL(NREG,NGRP)) +* +* 3) DO THE TRACKING OF THE EXACT GEOMETRY FOR *NEWT* OPTION. + IF( LTRK.NE.0 )THEN + NC= NTOTCL - 3 + IF( IPRT.GE.1 )THEN + WRITE(IOUT,'(1H )') + IF( NC.EQ.0 )THEN + WRITE(IOUT,'(/38H NOW, TRACKING GEOMETRY WITH NO CYLIND, + > 2HER/)') + ELSEIF( NC.EQ.1 )THEN + WRITE(IOUT,'(/38H NOW, TRACKING GEOMETRY WITH ONE CYLIN, + > 3HDER/)') + ELSE + WRITE(IOUT,'(/28H NOW, TRACKING GEOMETRY WITH,I4, + > 10H CYLINDERS/)') NC + ENDIF + ENDIF + ALLOCATE(ICUR(NTOTCL),INCR(NTOTCL)) + ALLOCATE(CONV(NTOTCL),TRKBEG(NTOTCL),TRKDIR(NTOTCL)) +* +* 3.0) WRITE FIRST RECORDS OF THE UNNORMALIZED TRACKING FILE + IF( LTRK.EQ.1 )THEN + LINMAX= 2*NVOL + 10 + ELSE + LINMAX= 8*NANGL*(2*NVOL + 8) + ENDIF + ISPEC = LTRK-1 + NALBG = 6 + ALLOCATE(NUMERO(LINMAX)) + ALLOCATE(LENGHT(LINMAX),ANGLES(3*MXANGL),DENSTY(MXANGL)) +* + NRMV=1 + CALL XL3TI3( IPRT, NANGLE, DENUSR, ISYMM, ANGLES, DENSTY, + > NTOTCL, NEXTGE, MAXR, REMESH, LINMAX, RCUTOF, + > NSUR, NVOL, INDEL, MINDIM, MAXDIM, ICORD, + > INCR, ICUR, TRKBEG, CONV, TRKDIR, LENGHT, + > NUMERO, NPIJ, NGRP, SIGTAL, SWVOID, NORE, + > NRMV, VOLTRK, KEYMRG,-NSOUT, NREG, NPSYS, + > DBLPIJ ) +* + CALL XL3NTR( IPRT, NDIM, ISPEC, NS, NV, NORE, VOLSUR, KEYMRG, + > MATALB, NANGL, VOLTRK, DENSTY ) +* + CALL XL3SIG( NGRP, NBMIX, XSSIGT, ALBOLD, NPSYS, NGRP, -NSOUT, + > NREG, MATMRG, VOLMRG(NSOUT+2), SIGTAL, SIGVOL, + > SWVOID, SWNZBC) +* + NRMV=0 + CALL XL3TI3( IPRT, NANGLE, DENUSR, ISYMM, ANGLES, DENSTY, + > NTOTCL, NEXTGE, MAXR, REMESH, LINMAX, RCUTOF, + > NSUR, NVOL, INDEL, MINDIM, MAXDIM, ICORD , + > INCR, ICUR, TRKBEG, CONV, TRKDIR, LENGHT, + > NUMERO, NPIJ, NGRP, SIGTAL, SWVOID, NORE, + > NRMV, VOLTRK, KEYMRG, -NSOUT, NREG , NPSYS, + > DBLPIJ ) +* + CALL QIJCMP(NREG,-NSOUT,NPIJ,NGRP,NCOR,VOLMRG,SIGTAL,DBLPIJ, + > NPSYS) +*---- +* RENORMALIZE ALL ISOTROPIC PROBS WITH VARIOUS OPTIONS +*---- + DO 2060 IGRP=1,NGRP + IF(NPSYS(IGRP).EQ.0) GO TO 2060 + IF( NRENOR.EQ.1 )THEN +* +* NORMALIZATION USING GELBARD SCHEME + CALL PIJRGL(IPRT,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP)) + ELSEIF( NRENOR.EQ.2 )THEN +* +* NORMALIZATION WORKING ON DIAGONAL COEFFICIENTS + CALL PIJRDG(NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP) ) + ELSEIF( NRENOR.EQ.3 )THEN +* +* NORMALIZATION WORKING ON WEIGHT FACTORS TO KEEP DIAG = 0.0 + CALL PIJRNL(IPRT,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP)) + ELSEIF( NRENOR .EQ. 4 )THEN ! ATTENTION +* +* NORMALIZATION WORKING ON WEIGHT FACTORS ADDITIVE (HELIOS) + CALL PIJRHL(IPRT,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP)) + ENDIF + IF( IPRT.GE.ICPALL )THEN + WRITE(IOUT,'(1H )') + WRITE(IOUT,'(35H COLLISION PROBABILITIES OUTPUT: , + > 35H *BEFORE* ALBEDO REDUCTION )') + CALL PIJWPR(0,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP), + > SIGVOL(1,IGRP),1) + ENDIF + 2060 CONTINUE +*---- +* ELIMINATION OF SURFACES FOR PIJ +*---- + IF( SWNZBC )THEN + ALLOCATE(PSST(NSOUT*NSOUT),PSVT(NSOUT*NREG)) + ALLOCATE(MATRT(NSOUT)) + CALL LCMLEN(IPTRK,'BC-REFL+TRAN',ILONG,ITYPE) + IF(ILONG.EQ.NSOUT) THEN + CALL LCMGET(IPTRK,'BC-REFL+TRAN',MATRT) + ELSE + DO 130 ISOUT=1,NSOUT + MATRT(ISOUT)=ISOUT + 130 CONTINUE + ENDIF + DO 2080 IGRP=1,NGRP + IF(NPSYS(IGRP).EQ.0) GO TO 2080 + CALL PIJABC(NREG,NSOUT,NPIJ,SIGTAL(1,IGRP),MATRT, + > DBLPIJ(1,IGRP),PSST,PSVT) + 2080 CONTINUE +* + DEALLOCATE(MATRT) + DEALLOCATE(PSVT,PSST) + ENDIF +* + ALLOCATE(FFACT(NREG)) + DO 2090 IGRP=1,NGRP + IF(NPSYS(IGRP).EQ.0) GO TO 2090 + IF( IPRT.GE.ICPEND )THEN + WRITE(IOUT,'(1H )') + WRITE(IOUT,'(35H COLLISION PROBABILITIES OUTPUT: , + > 35H *AFTER* ALBEDO REDUCTION )') + CALL PIJWPR(1,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP), + > SIGVOL(1,IGRP),1) + ENDIF +*---- +* CHARGE PIJ MATRIX IN THE DRAGON SYMMETRIZED FORMAT +*---- + DO 160 IIN=1,NREG + IF(SIGTAL(NSOUT+IIN+1,IGRP).EQ.0.0) THEN + FFACT(IIN)=1.0 + ELSE + FFACT(IIN)=1.0/SIGTAL(NSOUT+IIN+1,IGRP) + ENDIF + 160 CONTINUE + CALL PIJD2R(NREG,NSOUT,DBLPIJ(1,IGRP),FFACT,.FALSE.,NELPIJ, + > NPIJ,PIJ(1,IGRP)) + 2090 CONTINUE + DEALLOCATE(FFACT) +* + DEALLOCATE(DENSTY,ANGLES,LENGHT,NUMERO,TRKDIR,TRKBEG,CONV, + > INCR,ICUR) + ENDIF + DEALLOCATE(INDEL,ICORD,MAXDIM,MINDIM,REMESH,DBLPIJ,SIGTAL,SIGVOL, + > VOLSUR,VOLTRK,KEYMRG,MATALB) +*---- +* CHECK IF SCATTERING REDUCTION IS REQUIRED +*---- + ALLOCATE(PCSCT(NREG,2*NREG)) + DO 3000 IGRP=1,NGRP + IF(NPSYS(IGRP).EQ.0) GO TO 3000 + LSKIP=.TRUE. + DO 200 IBM=1,NBMIX + LSKIP=LSKIP.AND.(XSSIGW(IBM,1,IGRP).EQ.0.0) + 200 CONTINUE +*---- +* COMPUTE THE SCATTERING-REDUCED CP MATRICES +*---- + IF(.NOT.LSKIP) THEN + CALL PIJSMD(IPRT,NBMIX,NREG,MATMRG(NSOUT+2),VOLMRG(NSOUT+2), + > XSSIGW(0,1,IGRP),XSSIGT(0,IGRP),LEAKSW,PIJ(1,IGRP), + > PCSCT,1) + DO 220 I=1,NREG + FACT=VOLMRG(NSOUT+I+1) + DO 210 J=1,NREG + INDPIJ=INDPOS(I,J) + PIJ(INDPIJ,IGRP)=REAL(PCSCT(I,J))*FACT + 210 CONTINUE + 220 CONTINUE + ENDIF + 3000 CONTINUE + DEALLOCATE(PCSCT,VOLMRG,MATMRG) + RETURN +*---- +* FORMAT +*---- + 6000 FORMAT(' *** SPACE REQUIRED FOR CP MATRICES = ',I10,' K ***') + 6001 FORMAT(' *** CP MATRICES ALLOCATED ',10X,' ***') + 6002 FORMAT( + > ' -----------------'/' RECOMPUTED PARAMETERS '/ + > ' NANGL =',I10 ,' (NUMBER OF TRACKING ANGLES)'/ + > ' ISYMM =',I10 ,' (TRACKING SYMMETRY FACTOR)'/ + > ' CUTOFX =',F10.5,' (CUTOFF FOR TRACK LENGTH)'/ + > ' DENS =',F10.5,' (TRACK DENSITY)'/ + > ' PCORN =',F10.5,' (CORNER DUPLICATION DISTANCE)'/ + > ' -----------------'/) + END -- cgit v1.2.3