summaryrefslogtreecommitdiff
path: root/Dragon/src/MCGSCR.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/MCGSCR.f')
-rw-r--r--Dragon/src/MCGSCR.f338
1 files changed, 338 insertions, 0 deletions
diff --git a/Dragon/src/MCGSCR.f b/Dragon/src/MCGSCR.f
new file mode 100644
index 0000000..2287848
--- /dev/null
+++ b/Dragon/src/MCGSCR.f
@@ -0,0 +1,338 @@
+*DECK MCGSCR
+ SUBROUTINE MCGSCR(IPTRK,KPSYS,IPMACR,IPRINT,N1,NG,NGEFF,KPN,K,
+ 1 NREG,NANI,NFUNL,M,LPS,KEYFLX,KEYCUR,NZON,NGIND,
+ 2 NCONV,MXSCR,EPSSCR,REBAL,PHIOUT,PHIIN,V,NPJJM,
+ 3 KEYANI,IDIR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Acceleration of inner iteration (SCR method).
+*
+*Copyright:
+* Copyright (C) 2002 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. Le Tellier
+*
+*Parameters: input
+* IPTRK pointer to the tracking LCM object.
+* KPSYS pointer array for each group properties.
+* IPMACR pointer to the macrolib LCM object.
+* IPRINT print parameter (equal to zero for no print).
+* N1 number of unknowns per group of the corrective system.
+* NG number of groups.
+* NGEFF number of groups to process.
+* KPN total number of unknowns in vectors SUNKNO and FUNKNO.
+* K total number of volumes for which specific values
+* of the neutron flux and reactions rates are required.
+* NREG number of volumes.
+* NANI scattering anisotropy (=1 for isotropic scattering).
+* NFUNL number of moments of the flux (in 2D NFUNL=NANI*(NANI+1)/2).
+* M number of material mixtures.
+* LPS dimension of PSJ.
+* KEYFLX position of flux elements in FI vector.
+* KEYCUR position of current elements in FI vector.
+* NZON index-number of the mixture type assigned to each volume.
+* NGIND index of the groups to process.
+* NCONV logical array of convergence status for each group (.TRUE.
+* not converged).
+* MXSCR maximum number of iterations for rebalancing system.
+* EPSSCR convergence criterion for rebalancing system.
+* REBAL type of acceleration (.TRUE. rebalancing ; .FALSE. inner
+* iterations acceleration).
+* PHIIN initial guess (for this iteration) of zonal scalar flux.
+* V volumes.
+* NPJJM second dimension of PJJ.
+* KEYANI 'mode to l' index l=KEYANI(nu).
+* IDIR direction of fundamental current for TIBERE with MoC
+* =0,1,2,3.
+*
+*Parameters: input/output
+* PHIOUT zonal scalar flux.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK,KPSYS(NGEFF),IPMACR
+ INTEGER N1,NGEFF,NG,IPRINT,KPN,K,NREG,NANI,NFUNL,M,LPS,
+ 1 KEYFLX(NREG,NFUNL),KEYCUR(*),NZON(K),NGIND(NGEFF),MXSCR,NPJJM,
+ 2 KEYANI(NFUNL),IDIR
+ REAL EPSSCR,PHIIN(KPN,NGEFF),V(N1)
+ DOUBLE PRECISION PHIOUT(KPN,NGEFF)
+ LOGICAL NCONV(NGEFF),REBAL
+*----
+* LOCAL VARIABLES
+*----
+ TYPE(C_PTR) JPMACR,KPMACR,JPSYS
+ DOUBLE PRECISION TEMP
+ CHARACTER*12 NGTYP
+ CHARACTER*12 NAMPJJ,NAMPSJ
+ INTEGER, TARGET, SAVE, DIMENSION(1) :: IDUMMY
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: NGINDV,NJJ,IJJ,IPOS
+ REAL, ALLOCATABLE, DIMENSION(:) :: XSCAT,MATR
+ REAL, ALLOCATABLE, DIMENSION(:) :: PJJ,PSJ
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: SC
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: AR,PSI
+*
+ TYPE(C_PTR) PJJIND_PTR,IS_PTR,JS_PTR
+ INTEGER, POINTER, DIMENSION(:) :: IS,JS
+ INTEGER, POINTER, DIMENSION(:,:) :: PJJIND
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(NGINDV(NG),SC(0:M,NANI),AR(KPN,NGEFF,2),PSI(KPN,NGEFF,2),
+ 1 PJJ(NREG*NPJJM),PSJ(LPS))
+ PSI(:KPN,:NGEFF,:2)=0.0D0
+ AR(:KPN,:NGEFF,:2)=0.0D0
+ CALL LCMGPD(IPTRK,'PJJIND$MCCG',PJJIND_PTR)
+ CALL C_F_POINTER(PJJIND_PTR,PJJIND,(/ NPJJM,2 /))
+ IF(N1.GT.NREG) THEN
+* recover (IS,JS) arrays
+* IS: arrays for surfaces neighbors
+* JS: JS(IS(ISOUT)+1:IS(ISOUT+1)) give the neighboring regions to
+* surface ISOUT.
+ CALL LCMGPD(IPTRK,'IS$MCCG',IS_PTR)
+ CALL LCMGPD(IPTRK,'JS$MCCG',JS_PTR)
+ CALL C_F_POINTER(IS_PTR,IS,(/ N1-NREG+1 /))
+ CALL C_F_POINTER(JS_PTR,JS,(/ LPS /))
+ ELSE
+ IS=>IDUMMY
+ JS=>IDUMMY
+ ENDIF
+ IF(REBAL) THEN
+ JPMACR=LCMGID(IPMACR,'GROUP')
+ ALLOCATE(NJJ(0:M),IJJ(0:M),IPOS(0:M),XSCAT(0:M*NG))
+ ENDIF
+ IF(IDIR .EQ.0) THEN
+ NAMPJJ='PJJ$MCCG'
+ NAMPSJ='PSJ$MCCG'
+ ELSEIF(IDIR .EQ. 1) THEN
+ NAMPJJ='PJJX$MCCG'
+ NAMPSJ='PSJX$MCCG'
+ ELSEIF(IDIR .EQ. 2) THEN
+ NAMPJJ='PJJY$MCCG'
+ NAMPSJ='PSJY$MCCG'
+ ELSE
+ NAMPJJ='PJJZ$MCCG'
+ NAMPSJ='PSJZ$MCCG'
+ ENDIF
+*----
+* CONSTRUCT NGINDV (index to pass from "NGEFF format" to "NG format").
+*----
+ NGINDV(:NG)=0
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+ IG=NGIND(II)
+ NGINDV(IG)=II
+ ENDIF
+ ENDDO
+*---
+* COMPUTE RESIDUAL OF THE PREVIOUS FREE ITERATION FOR RHS
+*---
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+ IG=NGIND(II)
+ JPSYS=KPSYS(II)
+ CALL LCMGET(JPSYS,'DRAGON-S0XSC',SC(0,1))
+ IF(REBAL) THEN
+ KPMACR=LCMGIL(JPMACR,IG)
+ CALL LCMGET(KPMACR,'NJJS00',NJJ(1))
+ CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
+ CALL LCMGET(KPMACR,'IPOS00',IPOS(1))
+ CALL LCMGET(KPMACR,'SCAT00',XSCAT(1))
+ ENDIF
+ CALL MCGFCR(IPRINT,IG,II,NG,NGEFF,KPN,N1,NREG,NANI,NFUNL,
+ 1 M,.FALSE.,KEYFLX,KEYCUR,NZON,NGINDV,REBAL,PHIOUT,
+ 2 PHIIN,SC,KEYANI,NJJ,IJJ,IPOS,XSCAT,AR(1,II,1))
+ ENDIF
+ ENDDO
+*---
+* GAUSS SEIDEL ITERATIVE APPROACH TO SOLVE THE REBALANCING SYSTEM
+*---
+ IF(REBAL) THEN
+ NGTYP='GAUSS-SEIDEL'
+ NFIRST=NGEFF+1
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+ IG=NGIND(II)
+ KPMACR=LCMGIL(JPMACR,IG)
+ CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
+ DO IBM=1,M
+ IF(IJJ(IBM).GT.IG) THEN
+ NFIRST=II ! first thermal group index in NGEFF format
+ GOTO 5
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ ELSE
+ NGTYP=' ONE-GROUP'
+ NFIRST=1
+ ENDIF
+ 5 CONTINUE
+*
+ IF(NANI.GT.1) ALLOCATE(MATR(NFUNL*(NFUNL+1)*NREG))
+ DO ITSCR=1,MXSCR
+ DO 20 II=1,NGEFF
+ IF(NCONV(II)) THEN
+ IF((II.LT.NFIRST).AND.(ITSCR.GT.1)) GOTO 20
+ IG=NGIND(II)
+ JPSYS=KPSYS(II)
+ CALL LCMGET(JPSYS,'DRAGON-S0XSC',SC(0,1))
+ CALL LCMGET(JPSYS,NAMPJJ,PJJ)
+ IF(REBAL) THEN
+ KPMACR=LCMGIL(JPMACR,IG)
+ CALL LCMGET(KPMACR,'NJJS00',NJJ(1))
+ CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
+ CALL LCMGET(KPMACR,'IPOS00',IPOS(1))
+ CALL LCMGET(KPMACR,'SCAT00',XSCAT(1))
+ ENDIF
+ DO I=1,NREG
+ IBM=NZON(I)
+ DO INU=1,NFUNL
+ IND=KEYFLX(I,INU)
+ AR(IND,II,2)=AR(IND,II,1)
+ ENDDO
+ IF(REBAL) THEN
+* rebalancing option on : contribution from others groups.
+ IF(IBM.GT.0) THEN
+ IND=KEYFLX(I,1)
+ JG=IJJ(IBM)
+ DO 10 JND=1,NJJ(IBM)
+ IF(JG.NE.IG) THEN
+ JJ=NGINDV(JG)
+ IF(JJ.GT.0) THEN
+ AR(IND,II,2)=AR(IND,II,2)+
+ 1 XSCAT(IPOS(IBM)+JND-1)*PSI(I,JJ,1)
+ ENDIF
+ ENDIF
+ JG=JG-1
+ 10 CONTINUE
+ ENDIF
+ ENDIF
+ ENDDO
+ IF(NANI.EQ.1) THEN
+ DO I=1,NREG
+ IBM=NZON(I)
+ IND=KEYFLX(I,1)
+ PSI(IND,II,1)=AR(IND,II,2)
+ 1 *PJJ(I)/(1.0-SC(IBM,1)*PJJ(I))
+ ENDDO
+ ELSE
+ CALL MCGSCS(KPN,K,NREG,M,NANI,NFUNL,NPJJM,KEYFLX,KEYANI,
+ 1 PJJIND,NZON,SC(0,1),PJJ,AR(1,II,2),PSI(1,II,1),MATR)
+ ENDIF
+ ENDIF
+ 20 CONTINUE
+ IF(REBAL) THEN
+ ERRSCR=0.0
+ DO II=NFIRST,NGEFF
+ IF(NCONV(II)) THEN
+ ERR1=0.0
+ ERR2=0.0
+ DO I=1,NREG
+ DO INU=1,NFUNL
+ IND=KEYFLX(I,INU)
+ TEMP1=REAL(ABS(PSI(IND,II,1)-PSI(IND,II,2)))
+ TEMP2=REAL(ABS(PSI(IND,II,1)))
+ ERR1=MAX(ERR1,TEMP1)
+ ERR2=MAX(ERR2,TEMP2)
+ PSI(IND,II,2)=PSI(IND,II,1)
+ ENDDO
+ ENDDO
+ IF(ERR2.GT.0.0) ERRSCR=MAX(ERRSCR,ERR1/ERR2)
+ ENDIF
+ ENDDO
+ IF(ERRSCR.LT.EPSSCR) GO TO 30
+ ENDIF
+ ENDDO
+ 30 CONTINUE
+ IF(NANI.GT.1) DEALLOCATE(MATR)
+*
+ IF((REBAL).AND.(IPRINT.GT.0)) THEN
+ IF(NFIRST.GT.1) WRITE(6,100) NGIND(1),NGIND(NFIRST-1),NGTYP
+ IF((MXSCR.GT.1).AND.(NFIRST.LE.NGEFF)) THEN
+ WRITE(6,200) NGTYP,ERRSCR,(ITSCR-1)
+ ENDIF
+ ELSE
+ IF(IPRINT.GT.1) WRITE(6,100) NGIND(1),NGIND(NGEFF),NGTYP
+ ENDIF
+*----
+* PERFORM THE CORRECTION
+*----
+* Flux Correction
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+ DO I=1,NREG
+ DO INU=1,NFUNL
+ IND=KEYFLX(I,INU)
+ PHIOUT(IND,II)=PHIOUT(IND,II)+PSI(IND,II,1)
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+* Current Correction
+ IF(N1.GT.NREG) THEN
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+ IG=NGIND(II)
+ JPSYS=KPSYS(II)
+ CALL LCMGET(JPSYS,NAMPSJ,PSJ)
+ CALL LCMGET(JPSYS,'DRAGON-S0XSC',SC(0,1))
+ IF(REBAL) THEN
+ KPMACR=LCMGIL(JPMACR,IG)
+ CALL LCMGET(KPMACR,'NJJS00',NJJ(1))
+ CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
+ CALL LCMGET(KPMACR,'IPOS00',IPOS(1))
+ CALL LCMGET(KPMACR,'SCAT00',XSCAT(1))
+ ENDIF
+ DO I=NREG+1,N1
+ IIS=I-NREG
+ INDC=KEYCUR(IIS)
+ DO J=IS(IIS)+1,IS(IIS+1)
+ MCUI=JS(J)
+ IND=KEYFLX(MCUI,1)
+ IBM=NZON(MCUI)
+ IF(IBM.GT.0) THEN
+ TEMP=SC(IBM,1)*(PHIOUT(IND,II)-PHIIN(IND,II))
+ IF(REBAL) THEN
+ JG=IJJ(IBM)
+ DO 40 JND=1,NJJ(IBM)
+ IF(JG.NE.IG) THEN
+ JJ=NGINDV(JG)
+ IF(JJ.GT.0) THEN
+ TEMP=TEMP+XSCAT(IPOS(IBM)+JND-1)
+ 1 *(PHIOUT(IND,JJ)-PHIIN(IND,JJ))
+ ENDIF
+ ENDIF
+ JG=JG-1
+ 40 CONTINUE
+ ENDIF
+ PHIOUT(INDC,II)=PHIOUT(INDC,II)+PSJ(J)*TEMP/V(I)
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ ENDIF
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ IF(REBAL) DEALLOCATE(XSCAT,IPOS,IJJ,NJJ)
+ DEALLOCATE(PSJ,PJJ,PSI,AR,SC,NGINDV)
+ RETURN
+*
+ 100 FORMAT(10X,11HSCR: GROUPS,I4,3H TO,I4,2H: ,A12,7H SCHEME)
+ 200 FORMAT(10X,24HSCR: UP-SCATTE. GROUPS: ,A12,17H ITERATIONS: PRC:,
+ 1 E9.2,2H (,I4,12H ITERATIONS))
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 19:39:05 +0000