Initial commit from Polytechnique Montreal

1 files changed, 205 insertions, 0 deletions

diff --git a/Dragon/src/XDRH13.f b/Dragon/src/XDRH13.f
new file mode 100644
index 0000000..607dfae
--- /dev/null
+++ b/Dragon/src/XDRH13.f
@@ -0,0 +1,205 @@
+*DECK XDRH13
+      SUBROUTINE XDRH13 (IR1,NMILG,NG,NSMAX,IQUAD,FRTM,NS,IDIL,MIXGR,
+     1 RS,FRACT,SIGMA,SIGMS,P1I,P1DI,P1KI,SIGA1)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Calculation of the reduced collision probabilities for the She-Liu-Shi
+* double heterogeneity model (part 1).
+*
+*Copyright:
+* Copyright (C) 2019 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. Chambon
+*
+*Parameters: input
+* IR1     number of elementary mixtures in the domain.
+* NMILG   number of composite mixtures in the domain.
+* NG      number of different kind of micro structures. A kind of
+*         micro structure is characterized by the radius of its
+*         micro volumes. All the micro volumes of the same kind
+*         should own the same nuclear properties in a given macro
+*         volume.
+* NSMAX   maximum number of volumes (tubes or shells) in each kind of
+*         micro structure.
+* IQUAD   quadrature parameter for the treatment of the micro volumes.
+*         if IQUAD < 0, lines with regular interval are applied.
+* FRTM    minimum volume fraction of the grain in the representative 
+*         volume for She-Liu-Shi models.
+* NS      number of volumes in each kind of micro structure.
+* IDIL    elementary mixture indices in the diluent of the composite
+*         mixtures.
+* MIXGR   elementary mixture indices in the micro structures.
+* RS      radius of the micro volumes.
+* FRACT   volumic fractions of the micro volumes.
+*
+*Parameters: input/output
+* SIGMA   total macroscopic cross sections in each mixture of the
+*         composite geometry.
+* SIGMS   scattering macroscopic cross sections in each mixture of the
+*         composite geometry.
+*
+*Parameters: output
+* P1I     non collision probability in subvolume with 1 grain type.    
+* P1KI    escape probability from layer k in subvolume 
+*         with 1 grain type.
+* P1DI    escape probability from matrix in subvolume 
+*         with 1 grain type.
+* SIGA1   output cross sections.
+*
+*Reference:
+* D. She, Z. Liu, and L. Shi, An Equivalent Homogenization Method for
+* Treating the Stochastic Media, Nucl. Sci. Eng., 185,  351-360 (2018)
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IR1,NMILG,NG,NSMAX,IQUAD,NS(NG),IDIL(NMILG),
+     1 MIXGR(NSMAX,NG,NMILG)
+      REAL FRTM,RS(NSMAX+1,NG),FRACT(NG,IR1+NMILG),SIGMA(0:IR1+NMILG),
+     1 SIGMS(0:IR1+NMILG),P1I(NG,NMILG),P1DI(NG,NMILG),
+     2 P1KI(NSMAX,NG,NMILG),SIGA1(NG,NMILG)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER NR,IBM,J,K,N,M,MIL,NSMAX2
+      REAL DILF,DX,DXFACT,DRMIN,X,EP,EPI1,EPI2,LM,LGAR,SIGMA1,SIGMS1,
+     > FRT,RMAX,SIGT,XI,FRTT
+      DOUBLE PRECISION P1,P1D
+      REAL RGAR(NSMAX+2),LR(NSMAX+1)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: SIG
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: P1K
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      NSMAX2=NSMAX+NG
+      ALLOCATE(SIG(NSMAX2),P1K(NSMAX2))
+*----
+*  COMPUTE THE EQUIVALENT TOTAL AND SCATTERING CROSS SECTIONS 
+*  IN COMPOSITE REGIONS
+*----
+      DXFACT=100.0
+      IF(IQUAD.LT.0) DXFACT=REAL(-IQUAD)
+      DO 180 IBM=1,NMILG
+        MIL=IR1+IBM
+        SIGT=SIGMA(IDIL(IBM))
+        DILF=1.0
+        SIGMA(IR1+IBM)=0.0
+        SIGMS(IR1+IBM)=0.0
+        DO 10 J=1,NG
+          DILF=DILF-FRACT(J,MIL)
+   10   CONTINUE
+        DO 130 J=1,NG
+          FRT=FRACT(J,MIL)
+          IF(FRT.LE.0.00001) GO TO 130
+          NR=NS(J)
+          DRMIN=RS(NS(J)+1,J)
+          RGAR(1)=0.0
+          DO 15 K=2,NS(J)+1
+            DRMIN=MIN(DRMIN,RS(K,J)-RS(K-1,J))
+            RGAR(K)=RS(K,J)
+   15     CONTINUE
+          DO 20 K=1,NS(J)
+            P1K(K)=0.0D0
+            SIG(K)=SIGMA(MIXGR(K,J,IBM))
+   20     CONTINUE
+          FRTT=1.0-DILF
+*         FRT too small -> additional ring of matrix
+          IF((1.0-DILF).LT.FRTM) THEN
+            NR=NR+1
+            RGAR(NR+1)=RGAR(NR)*(FRTM/FRT)**(1.0/3.0)   
+            SIG(NR)=SIGT
+            FRTT=FRTM
+          ENDIF
+          RMAX=RGAR(NR+1)
+          LGAR=4.0/3.0*RMAX/FRTT
+          P1=0.0D0
+          P1D=0.0D0
+          DX=DRMIN/DXFACT
+          XI=-0.5
+          K=1
+* integral over radius to compute collision prob.
+   30     XI=XI+1.0
+          X=DX*XI
+          IF (X.GT.RGAR(K+1)) K=K+1
+          IF (K.GT.NR) GO TO 100
+*    Ref 1): Eq 13-17 
+*     compute segment lengths 
+          LM=LGAR/2.0
+          DO 40 N=1,NR
+            IF (N.LT.K) THEN
+              LR(N)=0.0D0
+            ELSEIF (N.EQ.K) THEN
+              LR(N)=(RGAR(N+1)**2.0 - X**2.0)**0.5
+            ELSE 
+              LR(N)=(RGAR(N+1)**2.0 - X**2.0)**0.5
+     1              -(RGAR(N)**2.0 - X**2.0)**0.5
+            ENDIF
+            LM=LM-LR(N)  
+   40     CONTINUE
+*     Ref 1): Eq 18-19
+          EP=2*SIGT*LM
+          DO 50 N=1,NR
+            EP=EP+2*LR(N)*SIG(N)
+   50     CONTINUE
+          P1=P1+X*DX*EXP(-EP)
+          DO 70 N=K,NR 
+            EPI1=SIGT*LM
+            EPI2=SIGT*LM
+            DO 60 M=1,NR
+              IF (M.LT.N) THEN
+                EPI2=EPI2+2*LR(M)*SIG(M)
+              ELSEIF (M.EQ.N) THEN
+                EPI2=EPI2+LR(M)*SIG(M)
+              ELSE 
+                EPI1=EPI1+LR(M)*SIG(M)
+                EPI2=EPI2+LR(M)*SIG(M)
+              ENDIF
+   60       CONTINUE
+* bug
+            IF(N.GT.NSMAX2) CALL XABORT('XDRH13: NSMAX OVERFLOW.')
+            P1K(N)=P1K(N)+X*DX*(EXP(-EPI1)+EXP(-EPI2))*
+     1              (1.0D0-EXP(-LR(N)*SIG(N)))
+   70     CONTINUE
+          GO TO 30
+  100     CONTINUE
+          P1=P1*2/RMAX**2.0
+          P1I(J,IBM)=REAL(P1)
+          P1D=1.0D0-P1
+          DO 110 K=1,NS(J)
+            P1K(K)=P1K(K)*2/RMAX**2.0
+            P1KI(K,J,IBM)=REAL(P1K(K))
+*      collision prob. conservation, Ref 1): Eq 4
+            P1D=P1D-P1K(K)
+  110      CONTINUE
+          P1DI(J,IBM)=REAL(P1D)
+*     Ref 1): Eq 5
+          SIGMA1=REAL(-LOG(P1)/LGAR)
+          SIGA1(J,IBM)=SIGMA1
+          SIGMS1=REAL(P1D/(1.0-P1)*SIGMA1/SIGT*SIGMS(IDIL(IBM)))
+          DO 120 K=1,NS(J)
+            SIGMS1=REAL(SIGMS1+P1K(K)/(1.0-P1)*SIGMA1/SIG(K)*
+     1               SIGMS(MIXGR(K,J,IBM)))
+  120      CONTINUE
+*     Ref 1): Eq 26
+          SIGMA(IR1+IBM)=SIGMA(IR1+IBM)+SIGMA1*FRT/(1-DILF)
+          SIGMS(IR1+IBM)=SIGMS(IR1+IBM)+SIGMS1*FRT/(1-DILF)
+  130     CONTINUE
+  180   CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(P1K,SIG)
+      RETURN
+      END