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/MUSJJ2.f90 | 372 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 372 insertions(+) create mode 100644 Dragon/src/MUSJJ2.f90 (limited to 'Dragon/src/MUSJJ2.f90') diff --git a/Dragon/src/MUSJJ2.f90 b/Dragon/src/MUSJJ2.f90 new file mode 100644 index 0000000..7570bf1 --- /dev/null +++ b/Dragon/src/MUSJJ2.f90 @@ -0,0 +1,372 @@ +! +!----------------------------------------------------------------------- +! +!Purpose: +! Compute the neutron flux and interface currents using the current +! iteration method for the multicell surfacic approximation. +! +!Copyright: +! Copyright (C) 2025 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): A. Hebert +! +!Parameters: input +! IPAS total number of regions. +! NMCEL total number of cells in the domain. +! NMERGE total number of merged cells for which specific values +! of the neutron flux and reactions rates are required. +! Many cells with different position in the domain can +! be merged before the neutron flux calculation if they +! own the same generating cell (NMERGE.le.NMCEL). +! NGEN total number of generating cells. A generating cell is +! defined by its material and dimensions, irrespective of +! its position in the domain (NGEN.le.NMERGE). +! IJAT total number of distinct out-currents. +! NPIJ size of cellwise scattering-reduced collision probability matrices. +! NPIS size of cellwise scattering-reduced escape probability matrices. +! NPSS size of cellwise scattering-reduced transmission probability matrices. +! EPSJ stopping criterion for flux-current iterations. +! NUNKNO total number of unknowns in vectors SUNKNO and FUNKNO. +! NMIX nmber of out-currents (dimension of arrays MIX and DVX). +! NIFR nmber of in-currents (dimension of arrays IFR and ALB). +! SUNKNO input source vector. +! IMPX print flag (equal to 0 for no print). +! NMC_NODE offset of the first volume in each generating cell. +! NMC_SURF offset of the first boundary surface in each generating cell. +! IFR index-number of in-currents. +! ALB transmission/albedo associated with each in-current. +! INUM index-number of the merged cell associated to each cell. +! MIX index-number of out-currents. +! DVX weight associated with each out-current. +! Note: IFR, ALB, MIX and DVX contains information to rebuild +! the geometrical 'A' matrix. +! IGEN index-number of the generating cell associated with each +! merged cell. +! IMAC global merge index assigned to each node in the surfacic +! geometry. +! PIJW cellwise scattering-reduced collision probability matrices. +! PISW cellwise scattering-reduced escape probability matrices. +! PSJW cellwise scattering-reduced collision probability matrices +! for incoming neutrons. +! PSSW cellwise scattering-reduced transmission probability matrices. +! +!Parameters: input/output +! FUNKNO unknown vector. +! +!----------------------------------------------------------------------- +! +SUBROUTINE MUSJJ2(IPAS,NMCEL,NMERGE,NGEN,IJAT,NPIJ,NPIS,NPSS,EPSJ,NUNKNO, & + & NMIX,NIFR,FUNKNO,SUNKNO,IMPX,NMC_NODE,NMC_SURF,IFR,ALB,INUM,MIX,DVX, & + & IGEN,IMAC,PIJW,PISW,PSJW,PSSW) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + !---- + ! SUBROUTINE ARGUMENTS + !---- + INTEGER IPAS,NMCEL,NMERGE,NGEN,IJAT,NPIJ,NPIS,NUNKNO,NMIX,NIFR,IMPX, & + & NMC_NODE(NGEN+1),NMC_SURF(NGEN+1),IFR(NIFR),INUM(NMCEL),MIX(NMIX), & + & IGEN(NMERGE),IMAC(IPAS) + REAL EPSJ,FUNKNO(NUNKNO),SUNKNO(NUNKNO),ALB(NIFR),DVX(NMIX),PIJW(NPIJ), & + & PISW(NPIS),PSJW(NPIS),PSSW(NPSS) + !---- + ! LOCAL VARIABLES + !---- + REAL PIJ,PIS + LOGICAL LOGTES + PARAMETER (MAXIT=400,LACCFC=2,ICL1=3,ICL2=3) + !---- + ! ALLOCATABLE ARRAYS + !---- + INTEGER, DIMENSION(:), POINTER :: INDNMC + DOUBLE PRECISION, DIMENSION(:), POINTER :: CIT0 + DOUBLE PRECISION, DIMENSION(:,:), POINTER :: CITR,AITR + DOUBLE PRECISION, DIMENSION(:), POINTER :: WCURR + !---- + ! SCRATCH STORAGE ALLOCATION + !---- + ALLOCATE(INDNMC(NMERGE)) + ALLOCATE(CITR(3,IJAT),CIT0(IJAT),AITR(2,IJAT)) + ALLOCATE(WCURR(IJAT)) + ! + KNMC=0 + DO JKK=1,NMERGE + JKG=IGEN(JKK) + J2=NMC_NODE(JKG+1)-NMC_NODE(JKG) + INDNMC(JKK)=KNMC + KNMC=KNMC+J2 + ENDDO + ! + DO I=1,IJAT + WCURR(I)=1.0D0 + CIT0(I)=0.0D0 + CITR(1,I)=FUNKNO(IPAS+I) + ENDDO + !---- + ! COMPUTE PSJW * Q(*) CONTRIBUTION + !---- + DO IKK=1,NMERGE + IKG=IGEN(IKK) + I2=NMC_NODE(IKG+1)-NMC_NODE(IKG) + I3=NMC_SURF(IKG+1)-NMC_SURF(IKG) + IT=0 + DO IK=1,IKK-1 + IT=IT+(NMC_SURF(IGEN(IK)+1)-NMC_SURF(IGEN(IK))) + ENDDO + IPSJ=0 + DO IK=1,IKG-1 + IPSJ=IPSJ+(NMC_NODE(IK+1)-NMC_NODE(IK))*(NMC_SURF(IK+1)-NMC_SURF(IK)) + ENDDO + KNMC=INDNMC(IKK) + DO I=1,I2 + DO IC=1,I3 + JCC=MIX(IT+IC) + PBJ=PSJW(IPSJ+(I-1)*I3+IC) + CIT0(JCC)=CIT0(JCC)+PBJ*DVX(IT+IC)*SUNKNO(IMAC(KNMC+I)) + ENDDO + ENDDO + ENDDO + !---- + ! COMPUTE NORMALIZATION VECTOR WCURR + !---- + DO ICEL=1,NMCEL + IKK=INUM(ICEL) + IKG=IGEN(IKK) + J3=NMC_SURF(IKG+1)-NMC_SURF(IKG) + IT=0 + DO IK=1,IKK-1 + IT=IT+(NMC_SURF(IGEN(IK)+1)-NMC_SURF(IGEN(IK))) + ENDDO + IS=0 + DO IK=1,ICEL-1 + IS=IS+(NMC_SURF(IGEN(INUM(IK))+1)-NMC_SURF(IGEN(INUM(IK)))) + ENDDO + IPSS=0 + DO IK=1,IKG-1 + IPSS=IPSS+(NMC_SURF(IK+1)-NMC_SURF(IK))**2 + ENDDO + DO JC=1,J3 + J1=IFR(IS+JC) + DO IC=1,J3 + PSS=PSSW(IPSS+(JC-1)*J3+IC) + WCURR(J1)=WCURR(J1)-PSS*ALB(IS+JC)*DVX(IT+IC) + ENDDO + ENDDO + ENDDO + ! + ISTART=1 + TEST=0.0D0 + ITER=0 + 10 ITER=ITER+1 + IF(ITER.GT.MAXIT) THEN + WRITE(6,'(/47H MUSJJ2: *** WARNING *** MAXIMUM NUMBER OF ITER, & + & 15HATIONS REACHED.)') + GO TO 190 + ENDIF + IT3=MOD(ITER,3)+1 + IT2=MOD(ITER-1,3)+1 + IT1=MOD(ITER-2,3)+1 + CITR(IT3,:IJAT)=CIT0(:IJAT) + !---- + ! COMPUTE PSSW * J(-) CONTRIBUTION + !---- + DO ICEL=1,NMCEL + IKK=INUM(ICEL) + IKG=IGEN(IKK) + J3=NMC_SURF(IKG+1)-NMC_SURF(IKG) + IT=0 + DO IK=1,IKK-1 + IT=IT+(NMC_SURF(IGEN(IK)+1)-NMC_SURF(IGEN(IK))) + ENDDO + IS=0 + DO IK=1,ICEL-1 + IS=IS+(NMC_SURF(IGEN(INUM(IK))+1)-NMC_SURF(IGEN(INUM(IK)))) + ENDDO + IPSS=0 + DO IK=1,IKG-1 + IPSS=IPSS+(NMC_SURF(IK+1)-NMC_SURF(IK))**2 + ENDDO + DO JC=1,J3 + J1=IFR(IS+JC) + DO IC=1,J3 + J2=MIX(IT+IC) + PSS=PSSW(IPSS+(JC-1)*J3+IC) + CITR(IT3,J2)=CITR(IT3,J2)+PSS*ALB(IS+JC)*DVX(IT+IC)*CITR(IT2,J1) + ENDDO + ENDDO + ENDDO + !---- + ! NORMALIZATION + !---- + S1=0.0D0 + S2=0.0D0 + DO I=1,IJAT + S1=S1+WCURR(I)*CITR(IT3,I) + S2=S2+CIT0(I) + ENDDO + ZNORM=S2/S1 + IF(ZNORM.LT.0.0D0) ZNORM=1.0D0 + CITR(IT3,:IJAT)=CITR(IT3,:IJAT)*ZNORM + !---- + ! ONE/TWO PARAMETER ACCELERATION + !---- + ALP=1.0D0 + BET=0.0D0 + LOGTES=(1+MOD(ITER-ISTART,ICL1+ICL2).GT.ICL1) + IF(LOGTES) THEN + AITR(1,:IJAT)=CITR(IT3,:IJAT)-CITR(IT2,:IJAT) + AITR(2,:IJAT)=CITR(IT2,:IJAT)-CITR(IT1,:IJAT) + DO ICEL=1,NMCEL + IKK=INUM(ICEL) + IKG=IGEN(IKK) + J3=NMC_SURF(IKG+1)-NMC_SURF(IKG) + IT=0 + DO IK=1,IKK-1 + IT=IT+(NMC_SURF(IGEN(IK)+1)-NMC_SURF(IGEN(IK))) + ENDDO + IS=0 + DO IK=1,ICEL-1 + IS=IS+(NMC_SURF(IGEN(INUM(IK))+1)-NMC_SURF(IGEN(INUM(IK)))) + ENDDO + IPSS=0 + DO IK=1,IKG-1 + IPSS=IPSS+(NMC_SURF(IK+1)-NMC_SURF(IK))**2 + ENDDO + DO JC=1,J3 + J1=IFR(IS+JC) + DO IC=1,J3 + J2=MIX(IT+IC) + PSS=PSSW(IPSS+(JC-1)*J3+IC)*ALB(IS+JC)*DVX(IT+IC) + AITR(1,J2)=AITR(1,J2)-PSS*(CITR(IT3,J1)-CITR(IT2,J1)) + AITR(2,J2)=AITR(2,J2)-PSS*(CITR(IT2,J1)-CITR(IT1,J1)) + ENDDO + ENDDO + ENDDO + IF((LACCFC.EQ.1).OR.(MOD(ITER-ISTART,ICL1+ICL2).EQ.ICL1)) THEN + S1=0.0D0 + S2=0.0D0 + DO I=1,IJAT + S1=S1+(CITR(IT3,I)-CITR(IT2,I))*AITR(1,I) + S2=S2+AITR(1,I)*AITR(1,I) + ENDDO + IF(S2.EQ.0.0D0) THEN + ISTART=ITER+1 + ELSE + ALP=S1/S2 + IF(ALP.LE.0.0D0) THEN + ISTART=ITER+1 + ALP=1.0D0 + ENDIF + ENDIF + DO I=1,IJAT + CITR(IT3,I)=CITR(IT2,I)+ALP*(CITR(IT3,I)-CITR(IT2,I)) + ENDDO + ELSE IF(LACCFC.EQ.2) THEN + S1=0.0D0 + S2=0.0D0 + S3=0.0D0 + S4=0.0D0 + S5=0.0D0 + DO I=1,IJAT + S1=S1+(CITR(IT3,I)-CITR(IT2,I))*AITR(1,I) + S2=S2+AITR(1,I)*AITR(1,I) + S3=S3+(CITR(IT3,I)-CITR(IT2,I))*AITR(2,I) + S4=S4+AITR(1,I)*AITR(2,I) + S5=S5+AITR(2,I)*AITR(2,I) + ENDDO + DET=S2*S5-S4*S4 + IF(DET.EQ.0.0D0) THEN + ISTART=ITER+1 + ELSE + ALP=(S5*S1-S4*S3)/DET + BET=(S2*S3-S4*S1)/DET + IF(ALP.LE.0.0D0) THEN + ISTART=ITER+1 + ALP=1.0D0 + BET=0.0D0 + ENDIF + ENDIF + DO I=1,IJAT + CITR(IT3,I)=CITR(IT2,I)+ALP*(CITR(IT3,I)-CITR(IT2,I))+ & + & BET*(CITR(IT2,I)-CITR(IT1,I)) + ENDDO + ENDIF + ENDIF + !---- + ! CHECK THE CONVERGENCE ERROR + !---- + ERR1=0.0D0 + ERR2=0.0D0 + DO I=1,IJAT + ERR1=MAX(ERR1,ABS(CITR(IT3,I)-CITR(IT2,I))) + ERR2=MAX(ERR2,ABS(CITR(IT3,I))) + ENDDO + IF(IMPX.GT.3) WRITE(6,'(30H MUSJJ2: CURRENT ITERATION NB.,I4, & + & 7H ERROR=,1P,E10.3,5H OVER,E10.3,15H NORMALIZATION=,E10.3, & + & 14H ACCELERATION=,2E11.3,1H.)') ITER,ERR1,ERR2,ZNORM,ALP,BET/ALP + IF(ITER.EQ.1) TEST=ERR1/ERR2 + IF((ITER.GT.20).AND.(ERR1/ERR2.GT.TEST)) THEN + WRITE(6,'(/50H MUSJJ2: *** WARNING *** CONVERGENCE DIFFICULTIES.)') + GO TO 190 + ENDIF + IF(LOGTES.OR.(ERR1.GT.EPSJ*ERR2)) GO TO 10 + IF(IMPX.GT.2) WRITE(6,'(40H MUSJJ2: CURRENT CONVERGENCE AT ITERATIO, & + & 5HN NB.,I4,7H ERROR=,1P,E10.3,5H OVER,E10.3,1H.)') ITER,ERR1,ERR2 + ! + 190 FUNKNO(:IPAS)=0.0 + DO I=1,IJAT + FUNKNO(IPAS+I)=REAL(CITR(IT3,I)) + ENDDO + !---- + ! COMPUTE PISW * J(-) CONTRIBUTION + !---- + DO ICEL=1,NMCEL + IKK=INUM(ICEL) + IKG=IGEN(IKK) + I2=NMC_NODE(IKG+1)-NMC_NODE(IKG) + I3=NMC_SURF(IKG+1)-NMC_SURF(IKG) + IS=0 + DO IK=1,ICEL-1 + IS=IS+(NMC_SURF(IGEN(INUM(IK))+1)-NMC_SURF(IGEN(INUM(IK)))) + ENDDO + IPIS=0 + DO IK=1,IKG-1 + IPIS=IPIS+(NMC_NODE(IK+1)-NMC_NODE(IK))*(NMC_SURF(IK+1)-NMC_SURF(IK)) + ENDDO + KNMC=INDNMC(IKK) + DO J=1,I2 + DO JC=1,I3 + J1=IFR(IS+JC) + PIS=PISW(IPIS+(JC-1)*I2+J) + FUNKNO(IMAC(KNMC+J))=FUNKNO(IMAC(KNMC+J))+PIS*ALB(IS+JC)*FUNKNO(IPAS+J1) + ENDDO + ENDDO + ENDDO + !---- + ! COMPUTE PIJW * Q(*) CONTRIBUTION + !---- + DO IKK=1,NMERGE + IKG=IGEN(IKK) + I2=NMC_NODE(IKG+1)-NMC_NODE(IKG) + IPIJ=0 + DO IK=1,IKG-1 + IPIJ=IPIJ+(NMC_NODE(IK+1)-NMC_NODE(IK))**2 + ENDDO + KNMC=INDNMC(IKK) + DO I=1,I2 + DO J=1,I2 + PIJ=PIJW(IPIJ+(I-1)*I2+J) + FUNKNO(IMAC(KNMC+J))=FUNKNO(IMAC(KNMC+J))+PIJ*SUNKNO(IMAC(KNMC+I)) + ENDDO + ENDDO + ENDDO +!---- +! SCRATCH STORAGE DEALLOCATION +!---- + DEALLOCATE(WCURR) + DEALLOCATE(AITR,CIT0,CITR) + DEALLOCATE(INDNMC) + RETURN + END -- cgit v1.2.3