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 --- Trivac/src/NSSFL5.f | 404 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 404 insertions(+) create mode 100755 Trivac/src/NSSFL5.f (limited to 'Trivac/src/NSSFL5.f') diff --git a/Trivac/src/NSSFL5.f b/Trivac/src/NSSFL5.f new file mode 100755 index 0000000..f6be4f3 --- /dev/null +++ b/Trivac/src/NSSFL5.f @@ -0,0 +1,404 @@ +*DECK NSSFL5 + SUBROUTINE NSSFL5(IPFLX,NUN,NG,NX,NY,NZ,LL4F,LL4X,LL4Y,LL4Z,NMIX, + > NALB,ICL1,ICL2,NADI,EPSNOD,MAXNOD,EPSTHR,MAXTHR,EPSOUT,MAXOUT, + > MAT,XX,YY,ZZ,XXX,YYY,ZZZ,IDL,VOL,KN,IQFR,QFR,DIFF,SIGR,CHI,SIGF, + > SCAT,BETA,FD,BNDTL,NPASS,MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,IMPX) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Flux calculation for the analytic nodal method in Cartesian 3D +* geometry using the nodal correction iteration strategy. +* +*Copyright: +* Copyright (C) 2022 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 +* IPFLX nodal flux. +* NUN number of unknowns per energy group. +* NG number of energy groups. +* NX number of nodes in the X direction. +* NY number of nodes in the Y direction. +* NZ number of nodes in the Z direction. +* LL4F number of nodal flux unknowns. +* LL4X number of nodal X-directed net currents unknowns. +* LL4Y number of nodal Y-directed net currents unknowns. +* LL4Z number of nodal Z-directed net currents unknowns. +* NMIX number of mixtures in the nodal calculation. +* NALB number of physical albedos. +* ICL1 number of free iterations in one cycle of the inverse power +* method (used for thermal iterations). +* ICL2 number of accelerated iterations in one cycle. +* NADI number of inner ADI iterations. +* EPSNOD nodal correction epsilon. +* MAXNOD maximum number of nodal correction iterations. +* EPSTHR thermal iteration epsilon. +* MAXTHR maximum number of thermal iterations. +* EPSOUT convergence epsilon for the power method. +* MAXOUT maximum number of iterations for the power method. +* MAT material mixtures. +* XX mesh spacings in the X direction. +* YY mesh spacings in the Y direction. +* ZZ mesh spacings in the Z direction. +* XXX Cartesian coordinates along the X axis. +* YYY Cartesian coordinates along the Y axis. +* ZZZ Cartesian coordinates along the Z axis. +* IDL position of averaged fluxes in unknown vector. +* VOL node volumes. +* KN node-ordered interface net current unknown list. +* IQFR boundary condition information. +* QFR albedo function information. +* DIFF diffusion coefficients +* SIGR removal cross sections. +* CHI fission spectra. +* SIGF nu times fission cross section. +* SCAT scattering cross section. +* BETA albedos. +* FD discontinuity factors. +* BNDTL set to 'flat' or 'quadratic'. +* NPASS number of transverse current iterations. +* MUX X-oriented compressed storage mode indices. +* MUY Y-oriented compressed storage mode indices. +* MUZ Z-oriented compressed storage mode indices. +* IMAX X-oriented position of each first non-zero column element. +* IMAY Y-oriented position of each first non-zero column element. +* IMAZ Z-oriented position of each first non-zero column element. +* IPY Y-oriented permutation matrices. +* IPZ Z-oriented permutation matrices. +* IMPX edition flag. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPFLX + INTEGER NUN,NG,NX,NY,LL4F,LL4X,LL4Y,NMIX,NALB,ICL1,ICL2, + 1 NADI,MAXNOD,MAXTHR,MAXOUT,IMPX,MAT(NX*NY*NZ),IDL(NX*NY*NZ), + 2 KN(6,NX,NY,NZ),IQFR(6,NX,NY,NZ),NPASS,MUX(LL4F),MUY(LL4F), + 3 MUZ(LL4F),IMAX(LL4F),IMAY(LL4F),IMAZ(LL4F),IPY(LL4F),IPZ(LL4F) + REAL EPSNOD,EPSTHR,EPSOUT,XX(NX*NY*NZ),YY(NX*NY*NZ),ZZ(NX*NY*NZ), + 1 XXX(NX+1),YYY(NY+1),ZZZ(NZ+1),VOL(NX*NY*NZ),QFR(6,NX*NY*NZ), + 2 DIFF(NMIX,NG),SIGR(NMIX,NG),CHI(NMIX,NG),SIGF(NMIX,NG), + 3 SCAT(NMIX,NG,NG),BETA(NALB,NG,NG),FD(NMIX,6,NG,NG) + CHARACTER(LEN=12) :: BNDTL +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) JPFLX + INTEGER, PARAMETER :: NDIM=3 + REAL :: COEF(6),KEFF,KEFF_OLD +*---- +* ALLOCATABLE ARRAYS +*---- + REAL, ALLOCATABLE, DIMENSION(:) :: EVECT + REAL, ALLOCATABLE, DIMENSION(:,:) :: A11X,A11Y,A11Z + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: QFR2,DRIFT + REAL, POINTER, DIMENSION(:,:) :: SAVG +* + ALB(X)=0.5*(1.0-X)/(1.0+X) +*---- +* SCRATCH STORAGE ALLOCATION +*---- + NEL=NX*NY*NZ + N=LL4F*NG + ALLOCATE(QFR2(6,NEL,NG),EVECT(N)) + ALLOCATE(DRIFT(6,NEL,NG),SAVG(NUN,NG)) +*---- +* ALBEDO PROCESSING +*---- + QFR2(:6,:NEL,:NG)=0.0 + DO IG=1,NG + DO IQW=1,6 + DO I=1,NX + DO J=1,NY + DO K=1,NZ + IEL=(K-1)*NX*NY+(J-1)*NX+I + IALB=IQFR(IQW,I,J,K) + IF(IALB > 0) THEN + IF(IALB.GT.NALB) CALL XABORT('NSSFL5: BETA OVERFLOW.') + QFR2(IQW,IEL,IG)=QFR(IQW,IEL)*ALB(BETA(IALB,IG,IG)) + ELSE + QFR2(IQW,IEL,IG)=QFR(IQW,IEL) + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO +*---- +* INITIALIZATIONS +*---- + KEFF_OLD=0.0 + KEFF=1.0 + CALL LCMLEN(IPFLX,'FLUX',ILONG,ITYLCM) + IF(ILONG == 0) THEN + JPFLX=LCMLID(IPFLX,'FLUX',NG) + SAVG(:NUN,:NG)=1.0 + ELSE + JPFLX=LCMGID(IPFLX,'FLUX') + DO IG=1,NG + CALL LCMLEL(JPFLX,IG,ILONG,ITYLCM) + IF(ILONG /= NUN) CALL XABORT('NSSFL5: INVALID FLUX.') + CALL LCMGDL(JPFLX,IG,SAVG(1,IG)) + ENDDO + CALL LCMGET(IPFLX,'K-EFFECTIVE',KEFF) + ENDIF + CALL LCMLEN(IPFLX,'DRIFT',ILONG,ITYLCM) + IF(ILONG == 0) THEN + JPFLX=LCMLID(IPFLX,'DRIFT',6*NEL) + DRIFT(:6,:NEL,:NG)=0.0 + ELSE + JPFLX=LCMGID(IPFLX,'DRIFT') + DO IG=1,NG + CALL LCMLEL(JPFLX,IG,ILONG,ITYLCM) + IF(ILONG /= 6*NEL) CALL XABORT('NSSFL5: INVALID DRIFT.') + CALL LCMGDL(JPFLX,IG,DRIFT(1,1,IG)) + ENDDO + ENDIF + DO IND1=1,LL4F + DO IG=1,NG + EVECT((IG-1)*LL4F+IND1)=SAVG(IND1,IG) + ENDDO + ENDDO +*---- +* NODAL CORRECTION LOOP +*---- + NMAX=IMAX(LL4F) + NMAY=IMAY(LL4F) + NMAZ=IMAZ(LL4F) + ALLOCATE(A11X(NMAX,NG),A11Y(NMAY,NG),A11Z(NMAZ,NG)) + JTER=0 + SAVG(:NUN,:NG)=0.0 + IOFY=7*LL4F+LL4X + IOFZ=7*LL4F+LL4X+LL4Y + DO WHILE((ABS(KEFF_OLD-KEFF) >= EPSNOD).OR.(JTER==0)) + JTER=JTER+1 + IF(IMPX > 0) THEN + WRITE(6,'(36H NSSFL5: Nodal correction iteration=,I5)') + > JTER + ENDIF + IF(JTER > MAXNOD) THEN + WRITE(6,'(/22H ACCURACY AT ITERATION,I4,2H =,1P,E12.5)') + > JTER,ABS(KEFF_OLD-KEFF) + CALL XABORT('NSSFL5: NODAL ITERATION FAILURE') + ENDIF + ! + ! set CMFD matrices + IOF=0 + DO IG=1,NG + CALL NSSMXYZ(LL4F,NDIM,NX,NY,NZ,NMIX,MAT,XX,YY,ZZ,IDL,VOL, + > IQFR,QFR2(1,1,IG),DIFF(1,IG),DRIFT(1,1,IG),SIGR(1,IG), + > MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,A11X(1,IG),A11Y(1,IG), + > A11Z(1,IG)) + ENDDO + ! + ! CMFD power iteration + DELTA=ABS(KEFF_OLD-KEFF) + KEFF_OLD=KEFF + CALL NSSEIG(NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG,MAT,IDL,VOL, + > MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,CHI,SIGF,SCAT,A11X,A11Y,A11Z, + > EPSTHR,MAXTHR,NADI,EPSOUT,MAXOUT,ICL1,ICL2,ITER,EVECT,KEFF,IMPX) + IF(IMPX > 0) WRITE(6,10) JTER,KEFF,ITER,DELTA + IF(IMPX > 2) THEN + WRITE(6,'(1X,A)') 'NSSFL5: EVECT=' + IOF=0 + DO IG=1,NG + WRITE(6,'(1X,1P,14E12.4)') EVECT(IOF+1:IOF+LL4F) + IOF=IOF+LL4F + ENDDO + ENDIF + ! + ! begin construct SAVG + IF(NUN /= IOFZ+LL4Z) CALL XABORT('NSSFL5: INVALID NUN.') + DO IEL=1,LL4F + DO IG=1,NG + SAVG(IEL,IG)=EVECT((IG-1)*LL4F+IEL) + ENDDO + ENDDO + ! + ! one- and two-node anm relations + CALL NSSANM3(NUN,NX,NY,NZ,LL4F,LL4X,LL4Y,NG,BNDTL,NPASS,NMIX, + > IDL,KN,IQFR,QFR2,MAT,XXX,YYY,ZZZ,KEFF,DIFF,SIGR,CHI,SIGF,SCAT, + > FD,SAVG) + ! + ! compute new drift coefficients + DRIFT(:6,:NEL,:NG)=0.0 + DO IG=1,NG + DO K=1,NZ + DO J=1,NY + DO I=1,NX + IEL=(K-1)*NX*NY+(J-1)*NX+I + IND1=IDL(IEL) + IF(IND1 == 0) CYCLE + KK1=IQFR(1,I,J,K) + KK2=IQFR(2,I,J,K) + JXM=KN(1,I,J,K) ; JXP=KN(2,I,J,K) + JYM=KN(3,I,J,K) ; JYP=KN(4,I,J,K) + JZM=KN(5,I,J,K) ; JZP=KN(6,I,J,K) + CALL NSSCO(NX,NY,NZ,NMIX,I,J,K,MAT,XX,YY,ZZ,DIFF(1,IG), + > IQFR(1,I,J,K),QFR2(1,IEL,IG),COEF) + IF((KK1 < 0) .AND. (KK2 < 0)) THEN + DRIFT(1,IEL,IG)=-(SAVG(7*LL4F+JXM,IG)+COEF(1)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + DRIFT(2,IEL,IG)=-(SAVG(7*LL4F+JXP,IG)-COEF(2)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + ELSE IF(KK1 < 0) THEN + DRIFT(1,IEL,IG)=-(SAVG(7*LL4F+JXM,IG)+COEF(1)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + IND3=IDL((K-1)*NX*NY+(J-1)*NX+I+1) + IF(IND3 /= 0) DRIFT(2,IEL,IG)=-(SAVG(7*LL4F+JXP,IG)+ + > COEF(2)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+ + > SAVG(IND1,IG)) + ELSE IF(KK2 < 0) THEN + IND2=IDL((K-1)*NX*NY+(J-1)*NX+I-1) + IF(IND2 /= 0) DRIFT(1,IEL,IG)=-(SAVG(7*LL4F+JXM,IG)+ + > COEF(1)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+ + > SAVG(IND2,IG)) + DRIFT(2,IEL,IG)=-(SAVG(7*LL4F+JXP,IG)-COEF(2)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + ELSE + IND2=IDL((K-1)*NX*NY+(J-1)*NX+I-1) + IND3=IDL((K-1)*NX*NY+(J-1)*NX+I+1) + IF(IND2 /= 0) DRIFT(1,IEL,IG)=-(SAVG(7*LL4F+JXM,IG)+ + > COEF(1)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+ + > SAVG(IND2,IG)) + IF(IND3 /= 0) DRIFT(2,IEL,IG)=-(SAVG(7*LL4F+JXP,IG)+ + > COEF(2)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+ + > SAVG(IND1,IG)) + ENDIF + KK3=IQFR(3,I,J,K) + KK4=IQFR(4,I,J,K) + IF((KK3 < 0).AND.(KK4 < 0)) THEN + DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+COEF(3)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)-COEF(4)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + ELSE IF(KK3 < 0) THEN + DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+COEF(3)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + IND3=IDL((K-1)*NX*NY+J*NX+I) + IF(IND3 /= 0) DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)+ + > COEF(4)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+ + > SAVG(IND1,IG)) + ELSE IF(KK4 < 0) THEN + IND2=IDL((K-1)*NX*NY+(J-2)*NX+I) + IF(IND2 /= 0) DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+ + > COEF(3)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+ + > SAVG(IND2,IG)) + DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)-COEF(4)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + ELSE + IND2=IDL((K-1)*NX*NY+(J-2)*NX+I) + IND3=IDL((K-1)*NX*NY+J*NX+I) + IF(IND2 /= 0) DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+ + > COEF(3)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+ + > SAVG(IND2,IG)) + IF(IND3 /= 0) DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)+ + > COEF(4)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+ + > SAVG(IND1,IG)) + ENDIF + KK5=IQFR(5,I,J,K) + KK6=IQFR(6,I,J,K) + IF((KK5 < 0).AND.(KK6 < 0)) THEN + DRIFT(5,IEL,IG)=-(SAVG(IOFZ+JZM,IG)+COEF(5)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + DRIFT(6,IEL,IG)=-(SAVG(IOFZ+JZP,IG)-COEF(6)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + ELSE IF(KK5 < 0) THEN + DRIFT(5,IEL,IG)=-(SAVG(IOFZ+JZM,IG)+COEF(5)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + IND3=IDL(K*NX*NY+(J-1)*NX+I) + IF(IND3 /= 0) DRIFT(6,IEL,IG)=-(SAVG(IOFZ+JZP,IG)+ + > COEF(6)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+ + > SAVG(IND1,IG)) + ELSE IF(KK6 < 0) THEN + IND2=IDL((K-2)*NX*NY+(J-1)*NX+I) + IF(IND2 /= 0) DRIFT(5,IEL,IG)=-(SAVG(IOFZ+JZM,IG)+ + > COEF(5)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+ + > SAVG(IND2,IG)) + DRIFT(6,IEL,IG)=-(SAVG(IOFZ+JZP,IG)-COEF(6)* + > SAVG(IND1,IG))/SAVG(IND1,IG) + ELSE + IND2=IDL((K-2)*NX*NY+(J-1)*NX+I) + IND3=IDL(K*NX*NY+(J-1)*NX+I) + IF(IND2 /= 0) DRIFT(5,IEL,IG)=-(SAVG(IOFZ+JZM,IG)+ + > COEF(5)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+ + > SAVG(IND2,IG)) + IF(IND3 /= 0) DRIFT(6,IEL,IG)=-(SAVG(IOFZ+JZP,IG)+ + > COEF(6)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+ + > SAVG(IND1,IG)) + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO + IF(IMPX > 5) THEN + DO IG=1,NG + WRITE(6,'(28H NSSFL5: DRIFT COEFFICIENTS(,I5,2H):)') IG + DO IEL=1,NX*NY*NZ + WRITE(6,'(1P,I7,6E12.4)') IEL,DRIFT(:6,IEL,IG) + ENDDO + ENDDO + ENDIF + ENDDO + DEALLOCATE(A11Z,A11Y,A11X) +*---- +* END OF NODAL CORRECTION LOOP +*---- + IF(IMPX.GT.0) WRITE(6,20) KEFF,JTER + IF(IMPX > 2) THEN + WRITE(6,'(/21H NSSFL5: UNKNOWNS----)') + DO IG=1,NG + WRITE(6,'(14H NSSFL5: SAVG(,I4,2H)=)') IG + WRITE(6,'(1P,12E12.4)') SAVG(:LL4F,IG) + WRITE(6,'(19H X-BOUNDARY FLUXES:)') + WRITE(6,'(1P,12E12.4)') SAVG(LL4F+1:2*LL4F,IG) + WRITE(6,'(1P,12E12.4)') SAVG(2*LL4F+1:3*LL4F,IG) + WRITE(6,'(19H Y-BOUNDARY FLUXES:)') + WRITE(6,'(1P,12E12.4)') SAVG(3*LL4F+1:4*LL4F,IG) + WRITE(6,'(1P,12E12.4)') SAVG(4*LL4F+1:5*LL4F,IG) + WRITE(6,'(19H Z-BOUNDARY FLUXES:)') + WRITE(6,'(1P,12E12.4)') SAVG(5*LL4F+1:6*LL4F,IG) + WRITE(6,'(1P,12E12.4)') SAVG(6*LL4F+1:7*LL4F,IG) + WRITE(6,'(12H X-CURRENTS:)') + WRITE(6,'(1P,12E12.4)') SAVG(7*LL4F+1:IOFY,IG) + WRITE(6,'(12H Y-CURRENTS:)') + WRITE(6,'(1P,12E12.4)') SAVG(IOFY+1:IOFZ,IG) + WRITE(6,'(12H Z-CURRENTS:)') + WRITE(6,'(1P,12E12.4)') SAVG(IOFZ+1:NUN,IG) + WRITE(6,'(5H ----)') + ENDDO + ENDIF +*---- +* SAVE SOLUTION +*---- + JPFLX=LCMGID(IPFLX,'FLUX') + DO IG=1,NG + CALL LCMPDL(JPFLX,IG,NUN,2,SAVG(1,IG)) + ENDDO + JPFLX=LCMGID(IPFLX,'DRIFT') + DO IG=1,NG + CALL LCMPDL(JPFLX,IG,6*NEL,2,DRIFT(1,1,IG)) + ENDDO + CALL LCMPUT(IPFLX,'K-EFFECTIVE',1,2,KEFF) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(EVECT,QFR2) + DEALLOCATE(SAVG,DRIFT) + RETURN +* + 10 FORMAT(14H NSSFL5: JTER=,I4,11H CMFD KEFF=,1P E13.6, + 1 12H OBTAINED IN,I4,28H CMFD ITERATIONS WITH ERROR=, + 2 1P,E11.4,1H.) + 20 FORMAT(18H NSSFL5: ANM KEFF=,F11.8,12H OBTAINED IN,I5, + 1 29H NODAL CORRECTION ITERATIONS.) + END -- cgit v1.2.3