summaryrefslogtreecommitdiff
path: root/Dragon/src/MCTRK.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/MCTRK.f')
-rw-r--r--Dragon/src/MCTRK.f399
1 files changed, 399 insertions, 0 deletions
diff --git a/Dragon/src/MCTRK.f b/Dragon/src/MCTRK.f
new file mode 100644
index 0000000..870468f
--- /dev/null
+++ b/Dragon/src/MCTRK.f
@@ -0,0 +1,399 @@
+*DECK MCTRK
+ SUBROUTINE MCTRK(IPTRK,IPRINT,NFREG,NFSUR,NDIM,NMIX,ANGBC,ITYPBC,
+ 1 MAXMSH,NUCELL,MXGSUR,MXGREG,MAXPIN,BCRT,ICODE,
+ 2 ALBEDO,IUNFLD,DGMESH,XYZL,INDEX,IDREG,DCMESH,
+ 3 ITPIN,DRAPIN,ISEED,NGRP,NL,NFM,NDEL,NED,MATCOD,
+ 4 XST,XSS,XSN2N,XSN3N,XSSNN,XSNUSI,XSCHI,XSEDI,
+ 5 MIX,ISONBR,NU,POS,ITALLY,NBSCO,NMERGE,NGCOND,
+ 6 IMERGE,INDGRP,SCORE1,SCORE2)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Simulation of a single particle from source to death.
+*
+*Copyright:
+* Copyright (C) 2008 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): B. Arsenault
+*
+*Parameters: input
+* IPTRK pointer to the TRACKING data structure.
+* IPRINT print level.
+* NFREG number of regions.
+* NFSUR number of surfaces.
+* NDIM problem dimensions.
+* NMIX number of mixtures in the geometry.
+* ANGBC angular treatment for boundary conditions (=0 isotropic;
+* =1 specular).
+* ITYPBC type of boundary.
+* MAXMSH maximum number of elements in MESH array.
+* NUCELL number of cell after unfolding in $X$, $Y$ and $Z$ directions.
+* MXGSUR maximum number of surfaces for any geometry.
+* MXGREG maximum number of region for any geometry.
+* MAXPIN maximum number of pins in a cell.
+* BCRT reflection/translation array.
+* ICODE albedo index array.
+* ALBEDO albedo array.
+* IUNFLD description of unfolded geometry.
+* DGMESH meshing vector for global geometry.
+* XYZL Cartesian boundary coordinates.
+* ISEED the seed for the generation of random numbers.
+* NGRP number of energy groups.
+* NL number of Legendre orders required in the estimations
+* (NL=1 or higher).
+* NFM number of fissile isotopes.
+* NDEL number of delayed precursor groups.
+* NED number of extra edit vectors.
+* XST total macroscopic cross sections for each mixture and energy
+* group.
+* XSS total scattering cross sections for each mixture and energy
+* group.
+* XSN2N N2N macroscopic cross sections for each mixture and energy
+* group.
+* XSN3N N3N macroscopic cross sections for each mixture and energy
+* group.
+* MATCOD region material.
+* XSSNN in-group and out-of-group macroscopic transfert cross sections
+* for each mixture.
+* XSNUSI the values of Nu time the fission cross sections for each
+* isotope per mixture and energy group.
+* XSCHI the values of fission spectrum per isotope per mixture for
+* each energy group.
+* XSEDI extra edit cross sections for each mixture and energy group.
+* MIX the mixture number where the fission occurs.
+* ISONBR the isotopic number where the fission occurs.
+* NU the value of the particle weight.
+* POS location of the particle in the x, y, and z directions.
+* ITALLY type of tally (=0 no tally; =1 score effective
+* multiplication factor; =2 also score macrolib information).
+* NBSCO number of macrolib related scores.
+* NMERGE number of homogenized regions.
+* NGCOND number of condensed energy groups.
+* IMERGE homogenized regions indices.
+* INDGRP condensed groups indices.
+* INDEX undefined.
+* IDREG undefined.
+* DCMESH undefined.
+* ITPIN undefined.
+* DRAPIN undefined.
+*
+*Parameters: input/output
+* SCORE1 score for total flux and effective multiplication factor.
+* SCORE2 macrolib score matrix.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+ IMPLICIT NONE
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK
+ INTEGER IPRINT,NFREG,NFSUR,NDIM,NMIX,ANGBC,ITYPBC,MAXMSH,
+ 1 NUCELL(3),MXGSUR,MXGREG,MAXPIN,MIX,ISONBR,
+ 2 BCRT(NFSUR),ICODE(6),IUNFLD(2,NUCELL(1),NUCELL(2),NUCELL(3)),
+ 3 INDEX(5,-MXGSUR:MXGREG,2),IDREG(MXGREG,2),ITPIN(3,MAXPIN),
+ 4 ISEED,NGRP,NL,NFM,NDEL,NED,MATCOD(NFREG),ITALLY,NBSCO,NMERGE,
+ 5 NGCOND,IMERGE(NFREG),INDGRP(NGRP)
+ REAL ALBEDO(6),XST(NMIX,NGRP),XSS(NMIX,NGRP,NL),XSN2N(NMIX,NGRP),
+ 1 XSN3N(NMIX,NGRP),XSSNN(NMIX,NGRP,NGRP,NL),
+ 2 XSNUSI(NMIX,NFM,NGRP,1+NDEL),XSCHI(NMIX,NFM,NGRP,1+NDEL),
+ 3 XSEDI(NMIX,NGRP,NED),NU,SCORE1(3),SCORE2(NBSCO,NMERGE,NGCOND)
+ DOUBLE PRECISION DGMESH(-1:MAXMSH,4),XYZL(2,NDIM),
+ 1 DCMESH(-1:MAXMSH,4,2),DRAPIN(-1:4,MAXPIN),POS(3)
+*----
+* LOCAL VARIABLES
+*----
+ INTEGER NBIND
+ PARAMETER (NBIND=7)
+ DOUBLE PRECISION XDRCST,PI,TWOPI
+ INTEGER INDX(NBIND,0:2),ODIR(3),IREG,IDIRC(2),MESHC(4,2),NSURC(2),
+ 1 NREGC(2),NTPIN,IDIR,ILEV,IDS,IDSO,IPSP,IEV,II,IGROUP,IGR,
+ 2 ISUR,ISUR2,IBM,IDIM,IFIRST
+ DOUBLE PRECISION CELLPO(3,2),PINCEN(3),LENGTH,MU,MUS,VIS,PHI,
+ 2 VDIR(3),OMEGA(3),ALB,NUSIGF0
+ LOGICAL TRACK,ACTIVE,KILLED,VIRTUAL
+ PARAMETER (ACTIVE=.TRUE.,KILLED=.FALSE.)
+ REAL RAND,XSM
+ DOUBLE PRECISION PSCAT,PN2N,PN3N,XPROB,XX
+*----
+* DATA VALUES
+*----
+ INTEGER INDOS(2,3),INDIC(-3:3)
+ DATA INDOS / 2,3,
+ 1 3,1,
+ 2 1,2 /
+ DATA INDIC / 5,3,1,0,2,4,6/
+*----
+* LOAD CONSTANTS
+*----
+ PI=XDRCST('Pi',' ')
+ TWOPI=2.D0*PI
+ IFIRST=1
+*----
+* TRACK THE NEUTRON PATH
+*----
+ DO ILEV=0,2
+ DO IDIR=1,3
+ INDX(IDIR,ILEV)=1
+ ENDDO
+ DO IDIR=4,NBIND
+ INDX(IDIR,ILEV)=0
+ ENDDO
+ ENDDO
+ ODIR(:3)=1
+
+ TRACK=ACTIVE
+ LENGTH=0.D0
+ IDSO=0
+ VIRTUAL=.FALSE.
+
+ DO WHILE(TRACK.EQV.ACTIVE)
+ IF(IDSO.EQ.0) THEN
+ IF(LENGTH.EQ.0.D0) THEN
+ IF(MIX.EQ.-1) THEN
+ IGR=1
+ ELSE
+ CALL RANDF(ISEED,IFIRST,RAND)
+ XPROB=0.D0
+ IGR=0
+ DO WHILE((RAND.GE.XPROB).AND.(IGR.LT.NGRP))
+ IGR=IGR+1
+ XPROB=XPROB+XSCHI(MIX,ISONBR,IGR,1)
+ ENDDO
+ ENDIF
+*----
+* SET THE TRACK DIRECTION FOR THE FIRST FLIGHT
+*----
+ XSM=0.0
+ DO II=1,NMIX
+ XSM=MAX(XSM,XST(II,IGR))
+ ENDDO
+ CALL RANDF(ISEED,IFIRST,RAND)
+ MU=2.D0*RAND-1.D0
+ MUS=DSQRT(1.D0-MU*MU)
+ CALL RANDF(ISEED,IFIRST,RAND)
+ PHI=TWOPI*RAND
+ VDIR(1)=MUS*COS(PHI)
+ VDIR(2)=MUS*SIN(PHI)
+ VDIR(3)=MU
+ ELSE
+*----
+* IF IT IS A VIRTUAL COLLISION, KEEP THE SAME DIRECTION,
+* IF IT IS NOT A VIRTUAL COLLISION IT CONSISTS OF AN ISOTROPIC
+* SCATTERING REACTION.
+*----
+ IF(.NOT.VIRTUAL) THEN
+ CALL RANDF(ISEED,IFIRST,RAND)
+ PHI=TWOPI*RAND
+ CALL RANDF(ISEED,IFIRST,RAND)
+ MU=COS(PI*(2.D0*RAND-1.D0))
+ MUS=DSQRT(1.D0-MU*MU)
+ VIS=DSQRT(1.D0-VDIR(1)*VDIR(1))
+ OMEGA(1)=VDIR(1)*MU-VIS*MUS*COS(PHI)
+ OMEGA(2)=VDIR(2)*MU+MUS*(VDIR(1)*VDIR(2)*
+ 1 COS(PHI)-VDIR(3)*SIN(PHI))/VIS
+ OMEGA(3)=VDIR(3)*MU+MUS*(VDIR(1)*VDIR(3)*
+ 1 COS(PHI)+VDIR(2)*SIN(PHI))/VIS
+ DO II=1,3
+ VDIR(II)=OMEGA(II)
+ ENDDO
+ ENDIF
+ ENDIF
+*----
+* SAMPLE THE FREE PATH DISTANCE
+*----
+ CALL RANDF(ISEED,IFIRST,RAND)
+ LENGTH=-LOG(RAND)/XSM
+ ENDIF
+*----
+* LOCATE THE NEUTRON IN THE GEOMETRY
+*----
+ CALL MCTPTR(IPTRK,IPRINT,NDIM,MAXMSH,ITYPBC,NUCELL,MXGSUR,
+ 1 MXGREG,MAXPIN,IUNFLD,DGMESH,XYZL,NBIND,POS,LENGTH,VDIR,
+ 2 ODIR,IDS,IDSO,IREG,INDX,IDIRC,MESHC,NSURC,NREGC,NTPIN,
+ 3 CELLPO,PINCEN,INDEX,IDREG,DCMESH,ITPIN,DRAPIN)
+*---
+* TALLY PROCESSING
+*---
+ IF(ITALLY.GT.0) THEN
+ CALL MCTALLY(ITALLY,NFREG,NMIX,NGRP,NL,NFM,NDEL,NED,NBSCO,
+ 1 NMERGE,NGCOND,IREG,IGR,NU,MATCOD,IMERGE,INDGRP,XSM,XST,XSS,
+ 2 XSN2N,XSN3N,XSSNN,XSNUSI,XSCHI,XSEDI,SCORE1,SCORE2)
+ ENDIF
+*----
+* AN INTERACTION HAS BEEN DETECTED, DETECT IF IT IS A VIRTUAL OR A
+* REAL COLLISION
+*----
+ IF(IREG.GT.0) THEN
+ IBM=MATCOD(IREG)
+ CALL RANDF(ISEED,IFIRST,RAND)
+ VIRTUAL=(RAND.LE.((XSM-XST(IBM,IGR))/XSM))
+*---
+* DETERMINE THE TYPE OF REACTION
+*---
+ IF(.NOT.VIRTUAL) THEN
+ CALL RANDF(ISEED,IFIRST,RAND)
+ PSCAT=XSS(IBM,IGR,1)/XST(IBM,IGR)
+ PN2N=XSN2N(IBM,IGR)/XST(IBM,IGR)
+ PN3N=XSN3N(IBM,IGR)/XST(IBM,IGR)
+ IF(RAND.LE.PSCAT+PN2N+PN3N) THEN
+*---
+* ISOTROPIC SCATTERING OR NxN EVENT
+*---
+ IF(RAND.GT.PSCAT+PN2N) THEN
+ NU=3.0*NU
+ ELSE IF(RAND.GT.PSCAT) THEN
+ NU=2.0*NU
+ ENDIF
+ XX = 0.D0
+ CALL RANDF(ISEED,IFIRST,RAND)
+ IGROUP=NGRP
+ DO WHILE((RAND.GT.XX).AND.(IGROUP.GE.1))
+ XX=XX+XSSNN(IBM,IGROUP,IGR,1)/XSS(IBM,IGR,1)
+ IGROUP=IGROUP-1
+ ENDDO
+ IGR=IGROUP+1
+ TRACK=ACTIVE
+ XSM=0.0
+ DO II=1,NMIX
+ XSM=MAX(XSM,XST(II,IGR))
+ ENDDO
+ ELSE
+*---
+* CAPTURE OR FISSION EVENT
+*---
+ TRACK = KILLED
+ MIX = IBM
+ NUSIGF0=0.0D0
+ DO II=1,NFM
+ NUSIGF0=NUSIGF0+XSNUSI(IBM,II,IGR,1)
+ ENDDO
+ IF((NUSIGF0.GT.0.0).AND.(NFM.EQ.1)) THEN
+ NU=NU*XSNUSI(IBM,1,IGR,1)/(XST(IBM,IGR)-XSS(IBM,IGR,1)-
+ 1 XSN2N(IBM,IGR)-XSN3N(IBM,IGR))
+ ISONBR=1
+ ELSE IF(NUSIGF0.GT.0.0) THEN
+ CALL RANDF(ISEED,IFIRST,RAND)
+ XPROB=0.D0
+ ISONBR=0
+ DO WHILE((RAND.GE.XPROB).AND.(ISONBR.LT.NFM))
+ ISONBR=ISONBR+1
+ XPROB=XPROB+XSNUSI(IBM,ISONBR,IGR,1)/NUSIGF0
+ ENDDO
+ NU=NU*XSNUSI(IBM,ISONBR,IGR,1)/(XST(IBM,IGR)-
+ 1 XSS(IBM,IGR,1)-2.0*XSN2N(IBM,IGR)-3.0*XSN3N(IBM,IGR))
+ ELSE
+ NU=0.0
+ ISONBR=0
+ ENDIF
+ ENDIF
+ ENDIF
+ ELSE
+*----
+* A BOUNDARY CONDITION HAS BEEN ENCOUNTERED
+*----
+ IF(ITYPBC.EQ.0) THEN
+*----
+* CARTESIAN BOUNDARY
+*----
+ ISUR=-IREG
+ ISUR2=BCRT(ISUR)
+ IREG=-ISUR2
+ IF((ISUR2.EQ.ISUR).AND.(ANGBC.EQ.1)) THEN
+* SPECULAR REFLECTIVE BOUNDARY CONDITION
+ IF(IPRINT.GT.4) WRITE(6,*) 'SPECULAR REFLECTION ON ',ISUR
+ ALB=ALBEDO(-ICODE(INDIC(IDSO)))
+ IF(ALB.EQ.0.0) THEN
+* no leakage
+ ISONBR = 0
+ TRACK = KILLED
+ ELSE IF(ALB.EQ.1.0) THEN
+ VDIR(IDS)=-VDIR(IDS)
+ ELSE
+ CALL RANDF(ISEED,IFIRST,RAND)
+ TRACK=(RAND.LE.ALB)
+ VDIR(IDS)=-VDIR(IDS)
+ ENDIF
+ ELSE IF(ISUR2.NE.ISUR) THEN
+* PERIODIC BOUNDARY CONDITION
+ IF(IPRINT.GT.4) WRITE(6,*) 'BC TRANSLATION FROM ',ISUR,
+ 1 ' TO',ISUR2
+ IPSP=-ISUR
+ IEV=-INDIC(IDSO)
+ IF(IPRINT.GT.99) CALL MCTPSP(IPTRK,POS,IPSP,IEV)
+ IF(IDSO.GT.0) THEN
+ POS(IDS)=POS(IDS)-(XYZL(2,IDS)-XYZL(1,IDS))
+ ELSE
+ POS(IDS)=POS(IDS)+(XYZL(2,IDS)-XYZL(1,IDS))
+ ENDIF
+ IDSO=-IDSO
+ ELSE IF((ISUR2.EQ.ISUR).AND.(ANGBC.EQ.0)) THEN
+* WHITE BOUNDARY CONDITION
+ ALB=ALBEDO(-ICODE(INDIC(IDSO)))
+ IF(ALB.EQ.0.0) THEN
+* no leakage
+ ISONBR = 0
+ TRACK = KILLED
+ ELSE
+* otherwise, choose randomly according to albedo
+ CALL RANDF(ISEED,IFIRST,RAND)
+ TRACK=(RAND.LE.ALB)
+ IF(TRACK.EQV.ACTIVE) THEN
+* FOR ISOTROPIC BC, CHOOSE RANDOMLY THE REENTERING
+* DIRECTION
+ CALL RANDF(ISEED,IFIRST,RAND)
+ MU=2.D0*RAND-1.D0
+ MUS=DSQRT(1.D0-MU*MU)
+ CALL RANDF(ISEED,IFIRST,RAND)
+ PHI=PI*RAND
+ VDIR(IDS)=-SIGN(1,IDSO)*MUS*SIN(PHI)
+ VDIR(INDOS(1,IDS))=MUS*COS(PHI)
+ VDIR(INDOS(2,IDS))=MU
+ DO IDIM=1,NDIM
+ IF(IDIM.NE.IDS) THEN
+ CALL RANDF(ISEED,IFIRST,RAND)
+ POS(IDIM)=XYZL(1,IDIM)+RAND*(XYZL(2,IDIM)-
+ 1 XYZL(1,IDIM))
+ ENDIF
+ ENDDO
+ CALL RANDF(ISEED,IFIRST,RAND)
+ LENGTH=-LOG(RAND)/XSM
+ ELSE
+ ISONBR = 0
+ ENDIF
+ ENDIF
+ ELSE
+ CALL XABORT('MCTRK: INVALID TYPE OF BOUNDARY CONDITION.')
+ ENDIF
+ ELSE
+*----
+* CYLINDRICAL BOUNDARY
+*----
+ CALL XABORT('MCTRK: CYLINDRICAL BOUNDARY NOT IMPLEMENTED.')
+ ENDIF
+ IF(IPRINT.GT.99) THEN
+ IPSP=-IREG
+ IEV=-INDIC(IDSO)
+ CALL MCTPSP(IPTRK,POS,IPSP,IEV)
+ ENDIF
+ ENDIF
+*----
+* SAVE NEUTRON PATHS IN NXT TABLE FOR PSP DISPLAY
+*----
+ IF(IPRINT.GT.99) THEN
+ IF(TRACK.EQV.KILLED) THEN
+ IPSP=-ABS(IREG)
+ ELSE
+ IPSP=ABS(IREG)
+ ENDIF
+ CALL MCTPSP(IPTRK,POS,IPSP,1)
+ ENDIF
+ ENDDO
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:32:44 +0000