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/BIVFIS.f | 186 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 186 insertions(+) create mode 100644 Dragon/src/BIVFIS.f (limited to 'Dragon/src/BIVFIS.f') diff --git a/Dragon/src/BIVFIS.f b/Dragon/src/BIVFIS.f new file mode 100644 index 0000000..49fe6da --- /dev/null +++ b/Dragon/src/BIVFIS.f @@ -0,0 +1,186 @@ +*DECK BIVFIS + SUBROUTINE BIVFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL, + > XSCHI,XSNUF,FUNKNO,SUNKNO) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute the fission source for a BIVAC tracking. +* +*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 +* IPTRK pointer to the tracking LCM object. +* NANIS maximum cross section Legendre order. +* NREG number of regions. +* NMAT number of mixtures. +* NIFIS number of fissile isotopes. +* NUNKNO number of unknowns per energy group. +* NGRP number of energy groups. +* MATCOD mixture indices. +* VOL volumes. +* XSCHI fission spectra. +* XSNUP nu times the fission cross sections. +* FUNKNO fluxes. +* +*Parameters: output +* SUNKNO sources. +* +*----------------------------------------------------------------------- +* + USE GANLIB + USE DOORS_MOD +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPTRK + INTEGER NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD(NREG) + REAL VOL(NREG),XSCHI(0:NMAT,NIFIS,NGRP),XSNUF(0:NMAT,NIFIS,NGRP), + 1 FUNKNO(NUNKNO,NGRP),SUNKNO(NUNKNO,NGRP) +*---- +* LOCAL VARIABLES +*---- + PARAMETER(NSTATE=40) + INTEGER JPAR(NSTATE),IJ1(25),IJ2(25) + LOGICAL CYLIND + CHARACTER*12 CXDOOR +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: KN,IDL + REAL, ALLOCATABLE, DIMENSION(:) :: XX,DD,FXSOR + REAL, ALLOCATABLE, DIMENSION(:,:) :: RR,RS +*---- +* RECOVER BIVAC SPECIFIC PARAMETERS. +*---- + CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR) + IF(JPAR(1).NE.NREG) CALL XABORT('BIVFIS: INCONSISTENT NREG.') + IF(JPAR(2).NE.NUNKNO) CALL XABORT('BIVFIS: INCONSISTENT NUNKNO.') + ITYPE=JPAR(6) + IELEM=JPAR(8) + ICOL=JPAR(9) + ISPLH=JPAR(10) + L4=JPAR(11) + LX=JPAR(12) + NLF=JPAR(14) + ISCAT=JPAR(16) + CYLIND=(ITYPE.EQ.3).OR.(ITYPE.EQ.6) + IF(ICOL.EQ.4) THEN + CALL XABORT('BIVFIS: COLLOCATION NODAL NOT IMPLEMENTED.') + ELSE IF((ITYPE.NE.2).AND.(ITYPE.NE.5)) THEN + CALL XABORT('BIVFIS: CARTESIAN 1D OR 2D GEOMETRY EXPECTED.') + ENDIF + CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM) + ALLOCATE(XX(NREG),DD(NREG),KN(MAXKN),IDL(NREG)) + CALL LCMGET(IPTRK,'XX',XX) + CALL LCMGET(IPTRK,'DD',DD) + CALL LCMGET(IPTRK,'KN',KN) + CALL LCMGET(IPTRK,'KEYFLX',IDL) + IF(IELEM.LT.0) THEN + ! Lagrangian finite element method +*---- +* RECOVER THE FINITE ELEMENT UNIT STIFFNESS MATRIX. +*---- + CALL LCMSIX(IPTRK,'BIVCOL',1) + CALL LCMLEN(IPTRK,'T',LC,ITYLCM) + ALLOCATE(RR(LC,LC),RS(LC,LC)) + CALL LCMGET(IPTRK,'R',RR) + CALL LCMGET(IPTRK,'RS',RS) + CALL LCMSIX(IPTRK,' ',2) +*---- +* COMPUTE VECTORS IJ1 AND IJ2 +*---- + LL=LC*LC + DO I=1,LL + IJ1(I)=1+MOD(I-1,LC) + IJ2(I)=1+(I-IJ1(I))/LC + ENDDO +*---- +* COMPUTE THE SOURCE +*---- + NUM1=0 + DO IR=1,NREG + IBM=MATCOD(IR) + IF(IBM.LE.0) CYCLE + IF(VOL(IR).EQ.0.0) GO TO 10 + DO I=1,LL + IND1=KN(NUM1+I) + IF(IND1.EQ.0) CYCLE + I1=IJ1(I) + I2=IJ2(I) + DO J=1,LL + IND2=KN(NUM1+J) + IF(IND2.EQ.0) CYCLE + IF(CYLIND) THEN + AUXX=(RR(I1,IJ1(J))+RS(I1,IJ1(J))*XX(IR)/DD(IR))* + > RR(I2,IJ2(J))*VOL(IR) + ELSE + AUXX=RR(I1,IJ1(J))*RR(I2,IJ2(J))*VOL(IR) + ENDIF + DO IG=1,NGRP + DO JG=1,NGRP + DO IS=1,NIFIS + SIGG=XSCHI(IBM,IS,IG)*XSNUF(IBM,IS,JG) + SUNKNO(IND1,IG)=SUNKNO(IND1,IG)+AUXX*SIGG* + > FUNKNO(IND2,JG) + ENDDO ! IS + ENDDO ! JG + ENDDO ! IG + ENDDO ! J + ENDDO ! I + 10 NUM1=NUM1+LL + ENDDO ! IR + DEALLOCATE(RS,RR) + ! append the integrated volumic sources + ALLOCATE(FXSOR(NUNKNO)) + DO IS=1,NIFIS + FXSOR(:NUNKNO)=0.0 + DO IR=1,NREG + IBM=MATCOD(IR) + IF(IBM.LE.0) CYCLE + DO IG=1,NGRP + FXSOR(IDL(IR))=FXSOR(IDL(IR))+VOL(IR)*XSNUF(IBM,IS,IG)* + > FUNKNO(IDL(IR),IG) + ENDDO ! IG + DO IG=1,NGRP + SUNKNO(IDL(IR),IG)=SUNKNO(IDL(IR),IG)+XSCHI(IBM,IS,IG)* + > FXSOR(IDL(IR)) + ENDDO ! IG + ENDDO ! IR + ENDDO ! IS + DEALLOCATE(FXSOR) + ELSE + ! Raviart-Thomas finite element method + CXDOOR='BIVAC' + ALLOCATE(FXSOR(NUNKNO)) + DO IS=1,NIFIS + FXSOR(:NUNKNO)=0.0 + DO IG=1,NGRP + CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,XSNUF(0,IS,IG), + > FXSOR,FUNKNO(1,IG)) + ENDDO + DO IR=1,NREG + IBM=MATCOD(IR) + IF(IBM.EQ.0) CYCLE + DO IE=1,IELEM**2 + IND=IDL(IR)+IE-1 + IF(IND.EQ.0) CYCLE + DO IG=1,NGRP + SUNKNO(IND,IG)=SUNKNO(IND,IG)+XSCHI(IBM,IS,IG)* + > FXSOR(IND) + ENDDO + ENDDO + ENDDO + ENDDO ! IS + DEALLOCATE(FXSOR) + ENDIF + DEALLOCATE(IDL,KN,DD,XX) + RETURN + END -- cgit v1.2.3