Initial commit from Polytechnique Montreal

1 files changed, 162 insertions, 0 deletions

diff --git a/Dragon/src/DOORS_MOD.f90 b/Dragon/src/DOORS_MOD.f90
new file mode 100644
index 0000000..9538892
--- /dev/null
+++ b/Dragon/src/DOORS_MOD.f90
@@ -0,0 +1,162 @@
+MODULE DOORS_MOD
+  USE GANLIB
+CONTAINS
+  SUBROUTINE DOORS(CDOOR,IPTRK,NMAT,NANIS,NUN,SIGG,SUNKNO,FUNKNO)
+    !
+    !---------------------------------------------------------------------
+    !
+    !Purpose:
+    ! compute the product of a cross section times a flux unknow vector.
+    !
+    !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
+    ! CDOOR   name of the geometry/solution operator.
+    ! IPTRK   pointer to the tracking (L_TRACK signature).
+    ! NMAT    number of mixtures in the macrolib.
+    ! NANIS   number of Legendre components in the macrolib (=0: isotropic).
+    ! NUN     total number of unknowns in vectors SUNKNO and FUNKNO.
+    ! SIGG    cross section.
+    ! FUNKNO  optional unknown vector. If not present, a flat flux
+    !         approximation is assumed.
+    !
+    !Parameters: output
+    ! SUNKNO  source vector. Volumes are included with BIVAC and TRIVAC
+    !         trackings.
+    !
+    !---------------------------------------------------------------------
+    !
+    !----
+    !  SUBROUTINE ARGUMENTS
+    !----
+    CHARACTER(LEN=12), INTENT(IN) :: CDOOR
+    TYPE(C_PTR), INTENT(IN) :: IPTRK
+    INTEGER, INTENT(IN) :: NMAT,NANIS,NUN
+    REAL, DIMENSION(0:NMAT,NANIS+1), INTENT(IN) :: SIGG
+    REAL, DIMENSION(NUN), INTENT(IN), OPTIONAL :: FUNKNO
+    REAL, DIMENSION(NUN), INTENT(INOUT) :: SUNKNO
+    !----
+    !  LOCAL VARIABLES
+    !----
+    INTEGER, PARAMETER :: NSTATE=40
+    INTEGER, DIMENSION(NSTATE) :: ISTATE
+    !----
+    !  ALLOCATABLE ARRAYS
+    !----
+    INTEGER, ALLOCATABLE, DIMENSION(:) :: MATCOD
+    INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: KEYFLX
+    REAL, ALLOCATABLE, DIMENSION(:) :: VOL
+    !----
+    !  RECOVER TRACKING PARAMETERS
+    !  NFUNL: number of spherical harmonics components used to expand the
+    !         flux and the sources.
+    !  NANIS_TRK: number of components in the angular expansion of the flux
+    !----
+    CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+    NREG=ISTATE(1)
+    IF(ISTATE(2).GT.NUN) CALL XABORT('DOORS: WRONG NUN.')
+    IF(ISTATE(4).GT.NMAT) CALL XABORT('DOORS: WRONG NMAT.')
+    NDIM=0
+    NLIN=1
+    NFUNL=1
+    NANIS_TRK=1
+    IF(CDOOR.EQ.'MCCG') THEN
+      NANIS_TRK=ISTATE(6)
+      NDIM=ISTATE(16)
+      CALL LCMGET(IPTRK,'MCCG-STATE',ISTATE)
+      NFUNL=ISTATE(19)
+      NLIN=ISTATE(20)
+    ELSE IF(CDOOR.EQ.'SN') THEN
+      NFUNL=ISTATE(7)
+      NLIN=ISTATE(8)
+      NDIM=ISTATE(9)
+      NLIN=NLIN**NDIM
+      NLIN=NLIN*ISTATE(35)
+      NANIS_TRK=ISTATE(16)
+    ELSE IF(CDOOR.EQ.'BIVAC') THEN
+      NLIN=ABS(ISTATE(8)) ! order of finite elements
+      NFUNL=MAX(1,ISTATE(14))
+      NANIS_TRK=ABS(ISTATE(16))
+    ELSE IF(CDOOR.EQ.'TRIVAC') THEN
+      NLIN=ABS(ISTATE(9)) ! order of finite elements
+      NLFUNL=MAX(1,ISTATE(30))
+      NANIS_TRK=ABS(ISTATE(32))
+    ENDIF
+    ALLOCATE(MATCOD(NREG),VOL(NREG),KEYFLX(NREG,NLIN,NFUNL))
+    KEYFLX(:NREG,:NLIN,:NFUNL)=0
+    CALL LCMLEN(IPTRK,'MATCOD',ILNLCM,ITYLCM)
+    IF(ILNLCM.NE.NREG) CALL XABORT('DOORS: INCOMPATIBLE NUMBER OF REGIONS.')
+    CALL LCMGET(IPTRK,'MATCOD',MATCOD)
+    CALL LCMGET(IPTRK,'VOLUME',VOL)
+    IF((CDOOR.EQ.'MCCG').OR.(CDOOR.EQ.'SN')) THEN
+      CALL LCMGET(IPTRK,'KEYFLX$ANIS',KEYFLX)
+    ELSE
+      CALL LCMGET(IPTRK,'KEYFLX',KEYFLX)
+    ENDIF
+    !----
+    !  PERFORM SIGG*FUNKNO MULTIPLICATION
+    !----
+    IF(CDOOR.EQ.'SN') THEN
+      CALL DOORS_SN(IPTRK,NANIS,NREG,NMAT,NUN,MATCOD,SIGG,SUNKNO,FUNKNO)
+    ELSE IF(CDOOR.EQ.'BIVAC') THEN
+      CALL DOORS_BIV(IPTRK,NANIS,NREG,NMAT,NUN,MATCOD,VOL,SIGG,SUNKNO,FUNKNO)
+    ELSE IF(CDOOR.EQ.'TRIVAC') THEN
+      CALL DOORS_TRI(IPTRK,NANIS,NREG,NMAT,NUN,MATCOD,VOL,SIGG,SUNKNO,FUNKNO)
+    ELSE IF(PRESENT(FUNKNO)) THEN
+      ! general case
+      IF((NANIS.EQ.0).OR.(NFUNL.EQ.1).OR.(NANIS_TRK.EQ.1)) THEN
+        ! LAB isotropy or transport correction
+        DO IR=1,NREG
+          IBM=MATCOD(IR)
+          IF(IBM.LE.0) CYCLE
+          DO IE=1,NLIN
+            IND=KEYFLX(IR,IE,1)
+            SUNKNO(IND)=SUNKNO(IND)+SIGG(IBM,1)*FUNKNO(IND)
+          ENDDO
+        ENDDO ! IR
+      ELSE
+        ! spherical harmonics expansion of the flux and source
+        DO IR=1,NREG
+          IBM=MATCOD(IR)
+          IF(IBM.LE.0) CYCLE
+          DO IAL=0,MIN(NFUNL-1,NANIS,NANIS_TRK-1)
+            FACT=REAL(2*IAL+1)
+            DO IAM=0,MIN(NFUNL-1,NANIS,NANIS_TRK-1)
+              DO IE=1,NLIN
+                IND=0
+                IF(NDIM.EQ.3) THEN
+                  IND=KEYFLX(IR,IE,1+IAL*NANIS_TRK+IAM)
+                ELSE IF((NDIM.EQ.2).AND.(IAM.LE.IAL)) THEN
+                  IND=KEYFLX(IR,IE,1+IAL*(IAL+1)/2+IAM)
+                ELSE IF(IAM.EQ.IAL) THEN
+                  IND=KEYFLX(IR,IE,1+IAL)
+                ENDIF
+                IF(IND.EQ.0) THEN
+                  CYCLE
+                ELSE IF(IND.GT.NUN) THEN
+                  CALL XABORT('DOORS: NUN OVERFLOW.')
+                ENDIF
+                SUNKNO(IND)=SUNKNO(IND)+FACT*SIGG(IBM,IAL+1)*FUNKNO(IND)
+              ENDDO ! IE
+            ENDDO ! IAM
+          ENDDO ! IAL
+        ENDDO ! IR
+      ENDIF
+    ELSE
+      ! general case (flat flux)
+      DO IR=1,NREG
+        IND=KEYFLX(IR,1,1)
+        SUNKNO(IND)=SUNKNO(IND)+SIGG(MATCOD(IR),1)
+      ENDDO
+    ENDIF
+    DEALLOCATE(KEYFLX,VOL,MATCOD)
+  END SUBROUTINE DOORS
+END MODULE DOORS_MOD
+