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/XL3TI3.f | 464 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 464 insertions(+) create mode 100644 Dragon/src/XL3TI3.f (limited to 'Dragon/src/XL3TI3.f') diff --git a/Dragon/src/XL3TI3.f b/Dragon/src/XL3TI3.f new file mode 100644 index 0000000..534b13c --- /dev/null +++ b/Dragon/src/XL3TI3.f @@ -0,0 +1,464 @@ +*DECK XL3TI3 + SUBROUTINE XL3TI3( IPRT,NANGLE,DENUSR,ISYMM,ANGLES,DENSTY, + > NTOTCL,NEXTGE,MAXR,REMESH,LINMAX,RCUTOF, + > NSUR,NVOL,INDEL,MINDIM, + > MAXDIM,ICOORD,INCR,ICUR,TRKBEG,CONV,TRKDIR, + > LENGHT,NUMERO,NPIJ,NGRP,SIGTAL,SWVOID,NORE, + > NRMV,VOLTRK,KEYMRG,NSOUT,NREG,NPSYS,DBLPIJ ) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Construct the sequential tape that will contain tracks for +* isotropic BC 3-D calculations. +* +*Copyright: +* Copyright (C) 1996 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 +* +*Parameters: input +* IPRT intermediate printing level for output. +* NANGLE number of angles used in the tracking process. +* DENUSR density of tracks in the plane perpendicular +* to the tracking angles. +* ISYMM flag for symetry (1/0 for on/off): +* = 2 reflection plane normal to X axis; +* = 4 reflection plane normal to Y axis; +* = 8 reflection plane normal to X and Y axis; +* =16 reflection plane normal to Z axis; +* =18 reflection plane normal to X and Z axis; +* =20 reflection plane normal to Y and Z axis; +* =24 reflection plane normal to X, Y and Z axis. +* ANGLES 3D angle values. +* DENSTY density of tracks angle by angle. +* NTOTCL number of cylindres of a type + 3. +* NEXTGE for tubez, nextge=1 +* MAXR max number of real mesh values in REMESH. +* REMESH real mesh values (rect/cyl). +* LINMAX max. number of track segments in a single track. +* RCUTOF cutof for corner tracking(0.25 suggested). +* NSUR number of surfaces. +* NVOL number of zones. +* INDEL numbering of surfaces and zones. +* MINDIM min index values for all axes (rect/cyl). +* MAXDIM max index values for all axes (rect/cyl). +* ICOORD principal axes direction (X/Y/Z) for meshes. +* ICUR current zonal location for a track segment. +* INCR increment direction for next track segment. +* TRKBEG position where a track begins. +* CONV s egments of tracks. +* TRKDIR direction of a track in all axes. +* LENGHT relative lenght of each segment in a track. +* NUMERO material identification of each track segment. +* NPIJ number of probabilities in one group +* NORE track normalization (-1 yes; 1 no). +* NRMV volume factors only (0 no; 1 yes). +* NGRP number of groups. +* SIGTAL total XS. +* SWVOID flag for void regions. +* KEYMRG merge keys. +* NSOUT number of outer surfaces. +* NREG number of regions. +* NPSYS undefined. +* DBLPIJ collision probabilities. +* +*Parameters: output +* VOLTRK volume factors. +* +*Reference: +* R.Roy, A. Hebert and G. Marleau +* A transport method for treating 3-d lattices of heterogeneous cells, +* Nuclear Science and Engineering, 101, 217 (1989). +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +* + INTEGER IPRT,NANGLE,NTOTCL,NEXTGE,MAXR,LINMAX, + > NSUR,NVOL,NPIJ,NGRP,NORE,NRMV,NSOUT,NREG + REAL DENUSR +* + REAL TRKBEG(NTOTCL), TRKDIR(NTOTCL), CONV(NTOTCL), + > REMESH(MAXR), DENSTY(*), RCUTOF, TRKEND(3), + > TRKORI(3), TRKOR2(3), PROJC2(3), ANGEQN(3,3), + > ANGLE2(3), ANGLE3(3), BARY(3), TRKCUT(3,2), + > TCUTOF(3,4), TORIC(3), ANGLES(3,*), + > SIGTAL(NSOUT:NREG,NGRP) + DOUBLE PRECISION LENGHT(LINMAX) + DOUBLE PRECISION VOLTRK(NSUR:NVOL,0:*),DBLPIJ(NPIJ,NGRP) + INTEGER MINDIM(NTOTCL), MAXDIM(NTOTCL), ICUR(NTOTCL), + > ICOORD(NTOTCL), INCR(NTOTCL), NUMERO(LINMAX), + > INDEL(4,*), NSBEG(4), NSEND(4), NSCUT(2), + > KEYMRG(NSUR:NVOL), NCSEG, NOLDS, NNEWS, NPSYS(NGRP) + INTEGER ISYMM + INTEGER IQUART(4) + LOGICAL LANGLE, SWVOID + EQUIVALENCE ( ANGEQN(1,2), ANGLE2 ), ( ANGEQN(1,3), ANGLE3 ) + CHARACTER TEDATA*13 + INTEGER NPAN, IOUT + DOUBLE PRECISION DZERO + PARAMETER ( NPAN=3, IOUT=6, DZERO=0.0D0 ) +* + INTEGER NDIM, NSOUTM, I, J, NPOINT, NPO2, NCUTOF, LINM2, + > NOTRAK, NSOLAN, IVS, IANG, IGRP, ISB, + > NANGLS, IPAN, NESTIM, JANG, IX, IY, IREF1, I2, I3, + > NSOLMX, NTTRK, IANG0, NDEBS, N, IZZ, IANGL, K, K3, + > NCROS, LINE, JC, IL + REAL ANORM2, A, B, DENLIN, RCIRC, R2CIRC, DP, SURTOT, + > VOLTOT, WEIGHT, DEPART, X, Y, ANN, ODDNXT, TOTLEN, + > TOTXXX, DDENST + DOUBLE PRECISION FCUTOF +* +* + ANORM2( A, B ) = A*A + B*B + NDIM= 3 + NSOUTM= -NSOUT + NSOLMX= 0 +* +* ONE WEIGHT FOR ALL TRACKS +* DENLIN= # OF TRACKS / CM + DENLIN= SQRT(DENUSR) +* +* COMPUTE THE CIRCUMSCRIBED RADIUS AND +* THE COORDINATES FOR THE TRUE CENTER OF THE CELL + R2CIRC= 0.0 + DO 10 I = 1, 3 + BARY(I)= 0.5 * (REMESH(MAXDIM(I)) + REMESH(MINDIM(I))) + IF( NEXTGE.EQ.1 )THEN + CALL XABORT('XL3TI3: TUBEZ NOT SUPPORTED') + ELSE + R2CIRC = R2CIRC + > + (REMESH(MAXDIM(I)) - REMESH(MINDIM(I)))**2 + ENDIF + 10 CONTINUE + R2CIRC= 0.25 * R2CIRC + RCIRC = SQRT(R2CIRC) +* +* NPOINT= # OF TRACKS ALONG ONE PERPENDICULAR AXIS + NPOINT= INT( 2. * RCIRC * DENLIN ) +***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** +***** CHANGE THIS "NPOINT" PARAMETER HAS TREMENDOUS EFFECTS ON TRACKING +***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** +* +* OTHER POSSIBLE CHOICES (EXPLORED WITHOUT SUCCESS) ARE ==> +*1-) NPOINT= INT( 2. * RCIRC * DENLIN ) + 1 +*2-) NPOINT= NINT( 2. * RCIRC * DENLIN ) +*3-) NPOINT= NINT( 2. * RCIRC * DENLIN ) + 1 +* +* KEEP "NPOINT" ODD & CORRECT DENSITY + NPO2 = NPOINT / 2 + NPOINT= 2 * NPO2 + 1 + DP = 2. * RCIRC / NPOINT + DENLIN= 1. / DP + DENUSR= DENLIN**2 + ODDNXT= (2*NPO2+3) / (2.*RCIRC) + IF( RCUTOF.EQ.0.0 )THEN + NCUTOF= 1 + ELSE + NCUTOF= 4 + ENDIF + LINM2= LINMAX-2*NCUTOF + FCUTOF= 1.0/DBLE(NCUTOF) + NOTRAK= 0 + SURTOT= 0.0 + VOLTOT= 0.0 + IQUART(1)=1 + IQUART(2)=1 + IQUART(3)=1 + IQUART(4)=1 + IF(NEXTGE .EQ. 1) THEN + IQUART(2)=0 + IQUART(3)=0 + IQUART(4)=0 + NSOLAN= (NANGLE * (NANGLE + 2)) / 2 + ELSE + IF( ISYMM .EQ. 8 .OR. ISYMM .EQ. 24 )THEN + NSOLAN= (NANGLE * (NANGLE + 2)) / 8 + IQUART(2)=0 + IQUART(3)=0 + IQUART(4)=0 + ELSE IF( ISYMM .EQ. 4 .OR. ISYMM .EQ. 20 )THEN + NSOLAN= (NANGLE * (NANGLE + 2)) / 4 + IQUART(2)=0 + IQUART(4)=0 + ELSE IF( ISYMM .EQ. 2 .OR. ISYMM .EQ. 18 )THEN + NSOLAN= (NANGLE * (NANGLE + 2)) / 4 + IQUART(3)=0 + IQUART(4)=0 + ELSE + NSOLAN= (NANGLE * (NANGLE + 2)) / 2 + ENDIF + ENDIF +* +* INITIALIZE NORMALIZED FACTORS + IF( NRMV.EQ.1 )THEN + DO 12 IVS= NSUR, NVOL + DO 11 IANG= 0, NSOLAN + VOLTRK(IVS,IANG)= DZERO + 11 CONTINUE + 12 CONTINUE + ELSE + DO 23 IGRP= 1, NGRP + DO 21 IVS= 1,NPIJ + DBLPIJ(IVS,IGRP)= DZERO + 21 CONTINUE + 23 CONTINUE + ENDIF + DDENST= 1.0/(NPAN*DENUSR) + NANGLS= (NANGLE * (NANGLE + 2)) / 2 + CALL XELEQN( 3, 0, ANGEQN ) + IANG= 0 + DO 15 JANG= 1, NANGLS + DO 16 IPAN= 1, NPAN + CALL XELEQN( 3, NANGLE, ANGEQN ) + 16 CONTINUE + IF(IQUART(MOD(JANG-1,4)+1).NE.1 ) GO TO 15 + IANG= IANG+1 + DENSTY(IANG)= REAL(2*NSOLAN) + ANGLES(1,IANG)= ANGEQN(1,1) + ANGLES(2,IANG)= ANGEQN(2,1) + ANGLES(3,IANG)= ANGEQN(3,1) + 15 CONTINUE +* +* TO REINITIATE THE EQN ANGLES + CALL XELEQN( 3, 0, ANGEQN ) + IF( NEXTGE.EQ.1 )THEN + DDENST= 12.0*DDENST + ENDIF + WEIGHT= 0.5*DDENST/REAL(NSOLAN) + IF( IPRT.GT.1 )THEN +* +* PREPARE & PRINT THE ESTIMATED NUMBER OF TRACKS + NESTIM= 0 + DEPART= - (NPO2+1) * DP + X = DEPART + DO 25 IX = 1, NPOINT + X = X + DP + Y = DEPART + DO 20 IY = 1, NPOINT + Y = Y + DP + IF( ANORM2( X, Y ) .LE. R2CIRC ) NESTIM= NESTIM + 1 + 20 CONTINUE + 25 CONTINUE + WRITE(IOUT,'(1H )') + WRITE(IOUT,'( 8H0ECHO = ,I8,20H TRACKS/AXIS/ANGLE )') + > NPOINT + WRITE(IOUT,'( 8H ECHO = ,I8,25H TRACKS/CIRCLE/ANGLE )') + > NESTIM + NESTIM= NESTIM * NPAN * NSOLAN + WRITE(IOUT,'(30H ECHO = NEXT ODD DENSITY > ,F15.7,4H/CM2)') + > ODDNXT**2 + WRITE(IOUT,'( 8H ECHO = ,28H ESTIMATED NUMBER OF TRACKS= ,I8 )') + > NESTIM +* +* PREPARE PRINTING WITH VARIABLE FORMAT + WRITE(IOUT,'(1H )') + WRITE(IOUT,'( 8H0ECHO = ,I3,27H SOLID ANGLES TO BE TRACKED )') + > NSOLAN + NSOLMX= MIN(9, NSOLAN/10) + IREF1 = 0 + WRITE(IOUT,'( 1H0,10(I1,9X))') (IREF1, IZZ=0,NSOLMX) + WRITE(IOUT,'( 1H ,10(I1,9X))') (MOD(IZZ,10), IZZ=0,NSOLMX) + WRITE(IOUT,'( 2H 0)') + TEDATA= '(1H+,TXXX,I1)' + ENDIF + IANG = 0 + IANG0 = 0 + NTTRK = 0 + DO 290 IANGL= 1, NANGLS + IANG=IANG+1 + LANGLE= .FALSE. + IF(IQUART(MOD(IANG-1,4)+1).NE.1)THEN +*---- +* Do not track this angle because of the problem symmetry +*---- + LANGLE= .FALSE. + ELSE +*---- +* Track this angle +*---- + IANG0= IANG0+1 + LANGLE=.TRUE. + ENDIF + IF(( IPRT.GT.1).AND.( MOD(IANGL,100) .EQ. 0 ))THEN + IREF1=IREF1+1 + NDEBS= NSOLMX+1 + NSOLMX=MIN(NDEBS+9, NANGLS/10) + WRITE(IOUT,'( 1H0,10(I1,9X))')(IREF1,IZZ=NDEBS,NSOLMX) + WRITE(IOUT,'( 1H ,10(I1,9X))') + > (MOD(IZZ,10),IZZ=NDEBS,NSOLMX) + WRITE(IOUT,'( 2H 0)') + ENDIF +* +* NPAN AXES DESCRIPTION (X=0.0, Y=0.0 & Z=0.0) + DO 250 IPAN= 1, NPAN + CALL XELEQN( 3, NANGLE, ANGEQN ) + IF(.NOT.LANGLE) GO TO 250 + IF( NEXTGE.EQ.1 )THEN + IF( IPAN.NE.2 ) GO TO 250 + ENDIF + DO 30 I = 1, 3 + N = ICOORD(I) + TRKDIR(N)= ANGEQN(N,1) + INCR(I)= +1 + IF( TRKDIR(N) .LT. 0.0 ) INCR(I)= -1 +* +* MODIFY ANGLES TO TAKE INTO ACCOUNT DP + ANGLE2(I)= DP * ANGLE2(I) + ANGLE3(I)= DP * ANGLE3(I) + IF( NCUTOF.NE.1 )THEN + TCUTOF(I,1)= RCUTOF*( ANGLE2(I)+ ANGLE3(I) ) + TCUTOF(I,2)= RCUTOF*( ANGLE2(I)- ANGLE3(I) ) + TCUTOF(I,3)= -TCUTOF(I,2) + TCUTOF(I,4)= -TCUTOF(I,1) + ENDIF +* +* DETERMINE THE ORIGIN OF ALL TRACKS + TRKOR2(I)= BARY(I) - (NPO2+1)*(ANGLE2(I)+ANGLE3(I)) + 30 CONTINUE + DO 45 I = 1, 3 + PROJC2(I)= 0.0 + DO 40 J = 1, 3 + IF( I.EQ.J ) GO TO 40 + PROJC2(I)= PROJC2(I) + TRKDIR(J) * TRKDIR(J) + 40 CONTINUE + 45 CONTINUE +* +* SCAN ALL TRACKS IN THE PERPENDICULAR PLANE + DO 180 I2 = 1, NPOINT + DO 50 J = 1, 3 + TRKOR2(J)= TRKOR2(J) + ANGLE2(J) + TRKORI(J)= TRKOR2(J) + 50 CONTINUE + DO 170 I3 = 1, NPOINT + ANN = 0.0 + DO 60 J = 1, 3 + TRKORI(J)= TRKORI(J) + ANGLE3(J) + ANN= ANN + (TRKORI(J)-BARY(J))**2 + 60 CONTINUE +* +* ELIMINATE TRACKS OUTSIDE CIRCUMSCRIBED CIRCLE + IF( ANN.GT.R2CIRC ) GO TO 170 +* +* WRITE(IOUT,7002) I2,I3,(TRKORI(JJ),JJ=1,3) +*7002 FORMAT(' ORIGINE MESH:',I10,5X,I10,5X,3(F11.5)) +* +* WHICH EXTERNAL SURFACES DO THIS TRACK CROSS ? + NTTRK=NTTRK+1 + CALL XELLSR( NDIM, NTOTCL, NSUR, MAXR, REMESH, + > INDEL, MINDIM, MAXDIM, ICOORD, ICUR, INCR, + > TRKORI, TRKDIR, TRKCUT, NSCUT, NCROS, + > TOTLEN) +* +* WHEN NOT SURFACES ARE CROSSED, ELIMINATE THE TRACK + IF(NCROS.LT.2) GO TO 170 + DO 70 K= 1, NDIM + TRKBEG(K)= TRKCUT(K,1) + TRKEND(K)= TRKCUT(K,2) + 70 CONTINUE + DO 75 K= 1, 4 + NSBEG(K)= NSCUT(1) + NSEND(K)= NSCUT(2) + 75 CONTINUE + IF( NCUTOF.NE.1 )THEN + DO 77 K= 1, 4 + DO 76 K3= 1, 3 + TORIC(K3)= TRKORI(K3)+TCUTOF(K3,K) + 76 CONTINUE + CALL XELLSR( NDIM, NTOTCL, NSUR, MAXR, REMESH, + > INDEL, MINDIM, MAXDIM, ICOORD, ICUR, INCR, + > TORIC, TRKDIR, TRKCUT, NSCUT, NCROS, + > TOTXXX) + IF(NSCUT(1).NE.0) NSBEG(K)= NSCUT(1) + IF(NSCUT(2).NE.0) NSEND(K)= NSCUT(2) + 77 CONTINUE + ENDIF + CALL XELLIN( NDIM, NTOTCL, MAXR, REMESH, + > NSUR, NVOL, INDEL, MINDIM, MAXDIM, + > ICOORD, ICUR, INCR, TRKBEG, TRKEND, TRKDIR, + > PROJC2, TOTLEN, CONV, LINM2, + > LENGHT(NCUTOF+1), NUMERO(NCUTOF+1), LINE) + NOTRAK= NOTRAK+1 +* + DO 78 ISB=1,NCUTOF + NUMERO(ISB)= NSBEG(ISB) + LENGHT(ISB)= FCUTOF + NUMERO(NCUTOF+LINE+ISB)= NSEND(ISB) + LENGHT(NCUTOF+LINE+ISB)= FCUTOF + 78 CONTINUE + LINE= LINE+2*NCUTOF + IF( IPRT.GT.10000)THEN + WRITE(IOUT,6001) NOTRAK, + > NCUTOF,(TRKBEG(JC),JC=1,3), + > NCUTOF,(TRKEND(JC),JC=1,3), + > (TRKDIR(JC),JC=1,3) + WRITE(IOUT,6002) (LENGHT(I),NUMERO(I),I=1,LINE) + ENDIF + IF( NRMV.EQ.1 )THEN + DO 301 I= 1, LINE + VOLTRK(NUMERO(I),IANG0)= VOLTRK(NUMERO(I),IANG0) + > + DBLE(WEIGHT) * LENGHT(I) + 301 CONTINUE + ELSE + IF( NORE.EQ.-1 )THEN + DO 302 I= 1, LINE + IF( NUMERO(I).GT.0 )THEN + LENGHT(I) = LENGHT(I) * + > SNGL(VOLTRK(NUMERO(I),IANG0)) + ENDIF + 302 CONTINUE + ENDIF + DO 303 I= 1, LINE + NUMERO(I)= KEYMRG(NUMERO(I)) + 303 CONTINUE +* +* START MODIFICATIONS 98/06 (G.M.) +* COMPRESS TRACKING FILE FOR +* SUCCESSIVE IDENTICAL REGIONS + NOLDS=NUMERO(1) + NCSEG=1 + DO 304 IL = 2, LINE + NNEWS=NUMERO(IL) + IF( NNEWS.LT.0 .OR. NNEWS.NE.NOLDS )THEN + NOLDS=NNEWS + NCSEG=NCSEG+1 + NUMERO(NCSEG)=NUMERO(IL) + LENGHT(NCSEG)=LENGHT(IL) + ELSEIF( NNEWS.EQ.NOLDS )THEN + LENGHT(NCSEG)=LENGHT(NCSEG)+LENGHT(IL) + ENDIF + 304 CONTINUE + CALL QIJI3D(NREG,NSOUTM,NPIJ,NGRP,LINMAX,NCUTOF, + > SWVOID,NCSEG,WEIGHT,NUMERO,LENGHT, + > SIGTAL,NPSYS,DBLPIJ) +* END MODIFICATIONS 98/06 (G.M.) +* + ENDIF + 170 CONTINUE + 180 CONTINUE + 250 CONTINUE + 290 CONTINUE +* + IF( IPRT.GT.1 )THEN + WRITE(IOUT,'(1H )') + WRITE(IOUT,'(27H0ECHO = TRACKING PROPERTIES )') + WRITE(IOUT,'( 8H0ECHO = ,I3,20H ANGLES AND DENSITY:, + > F9.6,4H/CM2)') + > NANGLE, DENUSR + WRITE(IOUT,'( 8H0ECHO = ,I10,3H / ,I10, + > 27H TRACKS NOT STORED ON TAPE /)') + > NOTRAK,NTTRK + ENDIF +* + RETURN +* + 6001 FORMAT(' #',I8,1P,' B',I1,'(',2(E10.2,','),E10.2,')', + > ' E',I1,'(',2(E10.2,','),E10.2,')', + > ' D(',2(E10.2,','),E10.2,')' ) + 6002 FORMAT(1P,5(1X,E15.7,1X,I6)) + END -- cgit v1.2.3