diff options
Diffstat (limited to 'Dragon/src/SNSBFP.f')
| -rw-r--r-- | Dragon/src/SNSBFP.f | 190 |
1 files changed, 190 insertions, 0 deletions
diff --git a/Dragon/src/SNSBFP.f b/Dragon/src/SNSBFP.f new file mode 100644 index 0000000..111cab0 --- /dev/null +++ b/Dragon/src/SNSBFP.f @@ -0,0 +1,190 @@ +*DECK SNSBFP + SUBROUTINE SNSBFP(IG,IPTRK,KPMACR,KPSYS,NANIS,NLF,NREG,NMAT, + 1 NUNKNO,NGRP,MATCOD,FLUX,QEXT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute the QEXT for the solution of SN equations with a Boltzmann- +* Fokker-Planck discretization. +* +*Copyright: +* Copyright (C) 2020 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 +* IG secondary group. +* IPTRK pointer to the tracking LCM object. +* KPMACR pointer to the secondary-group related macrolib information. +* KPSYS pointer to the system matrix information. +* NANIS maximum cross section Legendre order. +* NLF number of Legendre components in the flux. +* NREG number of regions. +* NMAT number of mixtures. +* NUNKNO number of unknowns per energy group including spherical +* harmonic terms, interface currents, fundamental currents +* and slowing-down angular fluxes at group boundary. +* NGRP number of energy groups. +* MATCOD mixture indices. +* FLUX fluxes and slowing-down angular fluxes at group boundary. +* +*Parameters: output +* QEXT sources and slowing-down angular fluxes at group boundary. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPTRK,KPMACR,KPSYS + INTEGER IG,NANIS,NLF,NREG,NMAT,NUNKNO,NGRP,MATCOD(NREG) + REAL FLUX(NUNKNO,NGRP),QEXT(NUNKNO,NGRP) +*---- +* LOCAL VARIABLES +*---- + PARAMETER(NSTATE=40,PI4=12.5663706144) + INTEGER JPAR(NSTATE),EELEM,P + CHARACTER CAN(0:19)*2 +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS + REAL, ALLOCATABLE, DIMENSION(:) :: XSCAT + TYPE(C_PTR) IL_PTR,IM_PTR + INTEGER, POINTER, DIMENSION(:) :: IL,IM +*---- +* DATA STATEMENTS +*---- + DATA CAN /'00','01','02','03','04','05','06','07','08','09', + > '10','11','12','13','14','15','16','17','18','19'/ +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(IJJ(0:NMAT),NJJ(0:NMAT),IPOS(0:NMAT)) + ALLOCATE(XSCAT(0:NMAT*NGRP)) +*---- +* RECOVER SNT SPECIFIC PARAMETERS +*---- + LFEP=0 + CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR) + IF(JPAR(1).NE.NREG) CALL XABORT('SNSBFP: INCONSISTENT NREG.') + IF(JPAR(2).NE.NUNKNO) CALL XABORT('SNSBFP: INCONSISTENT NUNKNO.') + IF(JPAR(15).NE.NLF) CALL XABORT('SNSBFP: INCONSISTENT NLF.') + ITYPE=JPAR(6) + NSCT=JPAR(7) + IELEM=JPAR(8) + ISCAT=JPAR(16) + EELEM=JPAR(35) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'IM',IM_PTR) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(IM_PTR,IM,(/ NSCT /)) +*---- +* CONSTRUCT THE QEXT. +*---- + IJJ(0)=0 + NJJ(0)=0 + IPOS(0)=0 + XSCAT(0)=0.0 + IOF0=0 + DO 130 P=1,NSCT + ILP=IL(P) + IF(ILP.GT.MIN(ISCAT-1,NANIS)) GO TO 130 + CALL LCMGET(KPMACR,'NJJS'//CAN(ILP),NJJ(1)) + CALL LCMGET(KPMACR,'IJJS'//CAN(ILP),IJJ(1)) + CALL LCMGET(KPMACR,'IPOS'//CAN(ILP),IPOS(1)) + CALL LCMGET(KPMACR,'SCAT'//CAN(ILP),XSCAT(1)) + IF((ITYPE.EQ.2).OR.(ITYPE.EQ.4)) THEN +*---- +* SLAB OR SPHERICAL 1D CASE. +*---- + NM=IELEM*EELEM + LFEP=IELEM*NLF*NREG + DO 20 IR=1,NREG + IBM=MATCOD(IR) + IF(IBM.LE.0) GO TO 20 + DO 15 IEL=1,NM + IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL + JG=IJJ(IBM) + DO 10 JND=1,NJJ(IBM) + IF(JG.NE.IG) THEN + QEXT(IND,IG)=QEXT(IND,IG)+FLUX(IND,JG)* + 1 XSCAT(IPOS(IBM)+JND-1) + ENDIF + JG=JG-1 + 10 CONTINUE + 15 CONTINUE + 20 CONTINUE + ELSE IF((ITYPE.EQ.5).OR.(ITYPE.EQ.6).OR.(ITYPE.EQ.8)) THEN +*---- +* 2D CASES (CARTESIAN OR R-Z). +*---- + NME=IELEM*IELEM + NM=NME*EELEM + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + LFEP=NME*NPQ*NREG + DO 70 IR=1,NREG + IBM=MATCOD(IR) + IF(IBM.LE.0) GO TO 70 + DO 65 IEL=1,NM + IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL + JG=IJJ(IBM) + DO 60 JND=1,NJJ(IBM) + IF(JG.NE.IG) THEN + QEXT(IND,IG)=QEXT(IND,IG)+FLUX(IND,JG)* + 1 XSCAT(IPOS(IBM)+JND-1) + ENDIF + JG=JG-1 + 60 CONTINUE + 65 CONTINUE + 70 CONTINUE + ELSE IF(ITYPE.EQ.7) THEN +*---- +* 3D CASES (CARTESIAN) +*---- + NME=IELEM*IELEM*IELEM + NM=NME*EELEM + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + LFEP=NME*NPQ*NREG + DO 110 IR=1,NREG + IBM=MATCOD(IR) + IF(IBM.LE.0) GO TO 110 + DO 100 IEL=1,NM + IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL + JG=IJJ(IBM) + DO 90 JND=1,NJJ(IBM) + IF(JG.NE.IG) THEN + QEXT(IND,IG)=QEXT(IND,IG)+FLUX(IND,JG)* + 1 XSCAT(IPOS(IBM)+JND-1) + ENDIF + JG=JG-1 + 90 CONTINUE + 100 CONTINUE + 110 CONTINUE + ELSE + CALL XABORT('SNSBFP: TYPE OF DISCRETIZATION NOT IMPLEMENTED.') + ENDIF + 130 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(XSCAT) + DEALLOCATE(IPOS,NJJ,IJJ) +*---- +* RECOVER SLOWING-DOWN ANGULAR FLUXES +*---- + CALL LCMGET(KPSYS,'DRAGON-DELTE',DELTAE) + IF(IG.EQ.1) THEN + QEXT(NUNKNO-LFEP+1:NUNKNO,1)=0.0 + ELSE + QEXT(NUNKNO-LFEP+1:NUNKNO,IG)=FLUX(NUNKNO-LFEP+1:NUNKNO,IG-1) + 1 *DELTAE + ENDIF + RETURN + END |