Initial commit from Polytechnique Montreal

1 files changed, 252 insertions, 0 deletions

diff --git a/Dragon/src/DOORS_TRI.f90 b/Dragon/src/DOORS_TRI.f90
new file mode 100644
index 0000000..4c17ebd
--- /dev/null
+++ b/Dragon/src/DOORS_TRI.f90
@@ -0,0 +1,252 @@
+SUBROUTINE DOORS_TRI(IPTRK,NANIS,NREG,NMAT,NUN,MATCOD,VOL,SIGG,SUNKNO,FUNKNO)
+  !
+  !-----------------------------------------------------------------------
+  !
+  !Purpose:
+  ! Compute the source for the solution of diffusion or PN equations.
+  ! Use a TRIVAC 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 (=0: isotropic).
+  ! NREG    number of regions.
+  ! NMAT    number of mixtures.
+  ! NUN     number of unknowns per energy group including net current.
+  ! MATCOD  mixture indices.
+  ! VOL     volumes. Volumes are included in SUNKNO.
+  ! SIGG    cross section.
+  ! FUNKNO  optional unknown vector. If not present, a flat flux
+  !         approximation is assumed.
+  !
+  !Parameters: input/output
+  ! SUNKNO  integrated sources.
+  !
+  !-----------------------------------------------------------------------
+  !
+  USE GANLIB
+  !----
+  !  SUBROUTINE ARGUMENTS
+  !----
+  TYPE(C_PTR) IPTRK
+  INTEGER NANIS,NREG,NMAT,NUN,MATCOD(NREG)
+  REAL VOL(NREG),SIGG(0:NMAT,NANIS+1),SUNKNO(NUN)
+  REAL, OPTIONAL :: FUNKNO(NUN)
+  !----
+  !  LOCAL VARIABLES
+  !----
+  PARAMETER(NSTATE=40)
+  INTEGER JPAR(NSTATE)
+  !----
+  !  RECOVER BIVAC SPECIFIC PARAMETERS.
+  !----
+  CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR)
+  IF(JPAR(1).NE.NREG) CALL XABORT('DOORS_TRI: INCONSISTENT NREG.')
+  IF(JPAR(2).NE.NUN) CALL XABORT('DOORS_TRI: INCONSISTENT NUN.')
+  IF(NANIS.NE.0) CALL XABORT('DOORS_TRI: SPN NOT IMPLEMENTED.')
+  ITYPE=JPAR(6)
+  IF((ITYPE.EQ.2).OR.(ITYPE.EQ.5).OR.(ITYPE.EQ.7)) THEN
+    ! Cartesian 1D, 2D or 3D geometry
+    CALL DOORS_TRIGSO(IPTRK,NREG,NMAT,NUN,MATCOD,VOL,SIGG,SUNKNO,FUNKNO)
+  ELSE IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9))  THEN
+    ! Hexagonal 2D or 3D geometry
+    CALL DOORS_TRIGSR(IPTRK,NREG,NMAT,NUN,MATCOD,SIGG,SUNKNO,FUNKNO)
+  ELSE
+    CALL XABORT('DOORS_TRI: GEOMETRY TYPE NOT AVAILABLE.')
+  ENDIF
+  RETURN
+CONTAINS
+  SUBROUTINE DOORS_TRIGSO(IPTRK,NREG,NMAT,NUN,MATCOD,VOL,SIGG,SUNKNO,FUNKNO)
+    !
+    !-----------------------------------------------------------------------
+    !
+    !Purpose:
+    ! Source term calculation for a mixed-dual formulation of the finite
+    ! element technique in a 3-D Cartesian geometry.
+    !
+    !-----------------------------------------------------------------------
+    !
+    USE GANLIB
+    !----
+    !  SUBROUTINE ARGUMENTS
+    !----
+    TYPE(C_PTR) IPTRK
+    INTEGER NREG,NMAT,NUN,MATCOD(NREG)
+    REAL VOL(NREG),SIGG(0:NMAT),SUNKNO(NUN)
+    REAL, OPTIONAL :: FUNKNO(NUN)
+    !----
+    !  LOCAL VARIABLES
+    !----
+    PARAMETER(NSTATE=40)
+    INTEGER IPAR(NSTATE)
+    CHARACTER HSMG*131
+    !----
+    !  ALLOCATABLE ARRAYS
+    !----
+    INTEGER, ALLOCATABLE, DIMENSION(:) :: KN
+    !----
+    !  RECOVER TRIVAC SPECIFIC PARAMETERS.
+    !----
+    CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR)
+    ITYPE=IPAR(6)
+    IELEM=IPAR(9)
+    ICOL=IPAR(10)
+    LX=IPAR(14)
+    LY=IPAR(15)
+    LZ=IPAR(16)
+    ISCAT=IPAR(32)
+    IF((ITYPE.NE.2).AND.(ITYPE.NE.5).AND.(ITYPE.NE.7)) THEN
+      CALL XABORT('DOORS_TRIGSO: INVALID CARTESIAN GEOMETRY.')
+    ELSE IF((IELEM.LT.0).OR.(ICOL.GT.3)) THEN
+      CALL XABORT('DOORS_TRIGSO: RAVIART-THOMAS METHOD EXPECTED.')
+    ELSE IF(ISCAT.GT.1) THEN
+      WRITE(HSMG,'(56HDOORS_TRIGSO: MACRO-CALCULATION WITH ANISOTROPIC SCATTER, &
+      & 66HING CURRENTLY NOT IMPLEMENTED; USE SCAT 1 KEYWORD IN TRIVAT: DATA.)')
+      CALL XABORT(HSMG)
+    ELSE IF(LX*LY*LZ.NE.NREG) THEN
+      CALL XABORT('DOORS_TRIGSO: INVALID NREG.')
+    ENDIF
+    CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
+    ALLOCATE(KN(MAXKN))
+    CALL LCMGET(IPTRK,'KN',KN)
+    IF(PRESENT(FUNKNO)) THEN
+      NUM1=0
+      DO K=1,NREG
+        L=MATCOD(K)
+        IF(L.LE.0) CYCLE
+        IF(VOL(K).EQ.0.0) GO TO 10
+        GARS=SIGG(L)
+        DO I0=1,IELEM**3
+          IND1=KN(NUM1+1)+I0-1
+          SUNKNO(IND1)=SUNKNO(IND1)+FUNKNO(IND1)*VOL(K)*GARS
+        ENDDO ! I0
+        10 NUM1=NUM1+1+6*IELEM**2
+      ENDDO ! K
+    ELSE
+      ! a flat flux is assumed
+      NUM1=0
+      DO K=1,NREG
+        L=MATCOD(K)
+        IF(L.LE.0) CYCLE
+        IF(VOL(K).EQ.0.0) GO TO 20
+        IND1=KN(NUM1+1)
+        SUNKNO(IND1)=SUNKNO(IND1)+VOL(K)*SIGG(L)
+        20 NUM1=NUM1+1+6*IELEM**2
+      ENDDO ! K
+    ENDIF
+    DEALLOCATE(KN)
+    RETURN
+  END SUBROUTINE DOORS_TRIGSO
+  !
+  SUBROUTINE DOORS_TRIGSR(IPTRK,NREG,NMAT,NUN,MATCOD,SIGG,SUNKNO,FUNKNO)
+    !
+    !-----------------------------------------------------------------------
+    !
+    !Purpose:
+    ! Source term calculation for a Thomas-Raviart-Schneider formulation of
+    ! the finite element technique in a 3-D hexagonal geometry.
+    !
+    !-----------------------------------------------------------------------
+    !
+    USE GANLIB
+    !----
+    !  SUBROUTINE ARGUMENTS
+    !----
+    TYPE(C_PTR) IPTRK
+    INTEGER NREG,NMAT,NUN,MATCOD(3,NREG/3)
+    REAL SIGG(0:NMAT),SUNKNO(NUN)
+    REAL, OPTIONAL :: FUNKNO(NUN)
+    !----
+    !  LOCAL VARIABLES
+    !----
+    PARAMETER(NSTATE=40)
+    INTEGER IPAR(NSTATE)
+    CHARACTER HSMG*131
+    !----
+    !  ALLOCATABLE ARRAYS
+    !----
+    INTEGER, ALLOCATABLE, DIMENSION(:) :: IPERT
+    INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KN
+    REAL, ALLOCATABLE, DIMENSION(:) :: FRZ
+    REAL, ALLOCATABLE, DIMENSION(:,:) :: ZZ
+    !----
+    !  RECOVER TRIVAC SPECIFIC PARAMETERS.
+    !----
+    CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR)
+    ITYPE=IPAR(6)
+    IELEM=IPAR(9)
+    ICOL=IPAR(10)
+    LX=IPAR(14)
+    LZ=IPAR(16)
+    ISCAT=IPAR(32)
+    NBLOS=LX*LZ/3
+    IF((ITYPE.NE.8).AND.(ITYPE.NE.9)) THEN
+      CALL XABORT('DOORS_TRIGSR: INVALID HEXAGONAL GEOMETRY.')
+    ELSE IF((IELEM.LT.0).OR.(ICOL.GT.3)) THEN
+      CALL XABORT('DOORS_TRIGSR: RAVIART-THOMAS METHOD EXPECTED.')
+    ELSE IF(ISCAT.GT.1) THEN
+      WRITE(HSMG,'(56HDOORS_TRIGSR: MACRO-CALCULATION WITH ANISOTROPIC SCATTER, &
+      & 66HING CURRENTLY NOT IMPLEMENTED; USE SCAT 1 KEYWORD IN TRIVAT: DATA.)')
+      CALL XABORT(HSMG)
+    ELSE IF(3*NBLOS.NE.NREG) THEN
+      CALL XABORT('DOORS_TRIGSR: INVALID NREG.')
+    ENDIF
+    CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
+    ALLOCATE(ZZ(3,NBLOS),KN(NBLOS,MAXKN/NBLOS),IPERT(NBLOS),FRZ(NBLOS))
+    CALL LCMGET(IPTRK,'SIDE',SIDE)
+    CALL LCMGET(IPTRK,'ZZ',ZZ)
+    CALL LCMGET(IPTRK,'KN',KN)
+    CALL LCMGET(IPTRK,'IPERT',IPERT)
+    CALL LCMGET(IPTRK,'FRZ',FRZ)
+    NELEM=IELEM*(IELEM+1)
+    TTTT=0.5*SQRT(3.0)*SIDE*SIDE
+    IF(PRESENT(FUNKNO)) THEN
+      NUM=0
+      DO KEL=1,NBLOS
+        IF(IPERT(KEL).EQ.0) CYCLE
+        NUM=NUM+1
+        L=MATCOD(1,IPERT(KEL))
+        IF(L.EQ.0) CYCLE
+        VOL0=TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL)
+        GARS=SIGG(L)
+        DO K3=0,IELEM-1
+          DO K2=0,IELEM-1
+            DO K1=0,IELEM-1
+              JND1=(NUM-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1
+              JND2=(KN(NUM,1)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1
+              JND3=(KN(NUM,2)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1
+              SUNKNO(JND1)=SUNKNO(JND1)+FUNKNO(JND1)*VOL0*GARS
+              SUNKNO(JND2)=SUNKNO(JND2)+FUNKNO(JND2)*VOL0*GARS
+              SUNKNO(JND3)=SUNKNO(JND3)+FUNKNO(JND3)*VOL0*GARS
+            ENDDO ! K1
+          ENDDO ! K2
+        ENDDO ! K3
+      ENDDO ! KEL
+    ELSE
+      ! a flat flux is assumed
+      NUM=0
+      DO KEL=1,NBLOS
+        IF(IPERT(KEL).EQ.0) CYCLE
+        NUM=NUM+1
+        L=MATCOD(1,IPERT(KEL))
+        IF(L.EQ.0) CYCLE
+        VOL0=TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL)
+        JND1=(NUM-1)*IELEM**3+1
+        JND2=(KN(NUM,1)-1)*IELEM**3+1
+        JND3=(KN(NUM,2)-1)*IELEM**3+1
+        SUNKNO(JND1)=SUNKNO(JND1)+VOL0*SIGG(L)
+        SUNKNO(JND2)=SUNKNO(JND2)+VOL0*SIGG(L)
+        SUNKNO(JND3)=SUNKNO(JND3)+VOL0*SIGG(L)
+      ENDDO ! KEL
+    ENDIF
+    DEALLOCATE(FRZ,IPERT,KN,ZZ)
+  END SUBROUTINE DOORS_TRIGSR
+END SUBROUTINE DOORS_TRI