Initial commit from Polytechnique Montreal

1 files changed, 267 insertions, 0 deletions

diff --git a/Dragon/src/BRERTS.f90 b/Dragon/src/BRERTS.f90
new file mode 100644
index 0000000..069a27b
--- /dev/null
+++ b/Dragon/src/BRERTS.f90
@@ -0,0 +1,267 @@
+SUBROUTINE BRERTS(IELEM,ICOL,NGRP,NLF,DELX,RCAT,JXM,JXP,IMPX,FHOMM,FHOMP)
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! Compute the Raviart-Thomas boundary fluxes for a single node in SPN
+! theory.
+!
+!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
+! IELEM   Raviart-Thomas polynomial order.
+! ICOL    Raviart-Thomas polynomial integration type.
+! NGRP    number of energy groups.
+! NLF     (NLF-1) is the SPN order (NLF is an even integer).
+! DELX    node width along X-axis.
+! RCAT    removal matrix (total minus scattering cross sections). The
+!         second dimension is for primary neutrons. The (2*IL+1) factor
+!         is included.
+! JXM     left boundary currents.
+! JXP     right boundary currents.
+! IMPX    print flag.
+!
+!Parameters: output
+! FHOMM   left boundary fluxes.
+! FHOMP   right boundary fluxes.
+!
+!-----------------------------------------------------------------------
+  USE GANLIB
+  !
+  !----
+  !  SUBROUTINE ARGUMENTS
+  !----
+  INTEGER, INTENT(IN) :: IELEM,ICOL,NGRP,NLF,IMPX
+  REAL, INTENT(IN) :: DELX
+  REAL, DIMENSION(NGRP,NLF/2), INTENT(IN) :: JXM,JXP
+  REAL(KIND=8), DIMENSION(NGRP,NGRP,NLF), INTENT(IN) :: RCAT
+  REAL, DIMENSION(NGRP,NLF/2), INTENT(OUT) :: FHOMM,FHOMP
+  !----
+  !  LOCAL VARIABLES
+  !----
+  TYPE(C_PTR) IPTRK
+  REAL QQ(5,5)
+  REAL(KIND=8) FHOMM_IG,FHOMP_IG
+  !----
+  !  ALLOCATABLE ARRAYS
+  !----
+  INTEGER, DIMENSION(:), ALLOCATABLE :: IND
+  REAL, DIMENSION(:,:), ALLOCATABLE :: R,V
+  REAL, DIMENSION(:,:,:), ALLOCATABLE :: FUNKNO
+  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: DXM,DXP,SYS
+  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: RCATI
+  !----
+  !  RECOVER FINITE ELEMENT UNIT MATRICES.
+  !----
+  CALL LCMOP(IPTRK,'***DUMMY***',0,1,0)
+  CALL BIVCOL(IPTRK,IMPX,IELEM,ICOL)
+  ALLOCATE(V(IELEM+1,IELEM),R(IELEM+1,IELEM+1))
+  CALL LCMSIX(IPTRK,'BIVCOL',1)
+  CALL LCMGET(IPTRK,'V',V)
+  CALL LCMGET(IPTRK,'R',R)
+  CALL LCMSIX(IPTRK,' ',2)
+  CALL LCMCL(IPTRK,2)
+  DO I0=1,IELEM
+    DO J0=1,IELEM
+      QQ(I0,J0)=0.0
+      DO K0=2,IELEM
+        QQ(I0,J0)=QQ(I0,J0)+V(K0,I0)*V(K0,J0)/R(K0,K0)
+      ENDDO
+    ENDDO
+  ENDDO
+  !----
+  !  INVERT THE RESIDUAL MATRIX.
+  !----
+  IF(MOD(NLF,2).NE.0) CALL XABORT('BRERTS: EVEN NLF EXPECTED.')
+  ALLOCATE(RCATI(NGRP,NGRP,NLF),IND(NGRP))
+  DO IL=2,NLF,2
+    RCATI(:NGRP,:NGRP,IL)=RCAT(:NGRP,:NGRP,IL)
+    CALL ALINVD(NGRP,RCATI(1,1,IL),NGRP,IER,IND)
+    IF(IER.NE.0) CALL XABORT('BRERTS: SINGULAR MATRIX(1).')
+  ENDDO
+  DEALLOCATE(IND)
+  !----
+  !  LEAKAGE-REMOVAL SYSTEM MATRIX ASSEMBLY FOR THE RAVIART-THOMAS METHOD.
+  !----
+  FHOMM(:NGRP,:NLF/2)=0.0
+  FHOMP(:NGRP,:NLF/2)=0.0
+  L4=IELEM*NGRP
+  INX=L4*NLF/2+1
+  ALLOCATE(SYS(L4*NLF/2,INX))
+  SYS(:L4*NLF/2,:INX)=0.0D0
+  DO IL=0,NLF-1
+    IF(MOD(IL,2).EQ.0) THEN
+      !----
+      !  EVEN PARITY EQUATION.
+      !----
+      DO IG=1,NGRP
+        DO JG=1,NGRP
+          DO J0=1,IELEM
+            JND1=(IL/2)*L4+(IG-1)*IELEM+J0
+            JND2=(IL/2)*L4+(JG-1)*IELEM+J0
+            SYS(JND1,JND2)=SYS(JND1,JND2)+DELX*RCAT(IG,JG,IL+1) ! IG <-- JG
+          ENDDO
+        ENDDO
+      ENDDO
+    ELSE
+      !----
+      !  ODD PARITY EQUATION.
+      !----
+      DO IG=1,NGRP
+        IF(IELEM.GT.1) THEN
+          ! get rid of net current collocation points inside finite elements.
+          DO JG=1,NGRP
+            DO J0=1,IELEM
+              JND1=((IL-1)/2)*L4+(IG-1)*IELEM+J0
+              DO K0=1,IELEM
+                IF(QQ(J0,K0).EQ.0.0) CYCLE
+                KND2=((IL-1)/2)*L4+(JG-1)*IELEM+K0
+                SYS(JND1,KND2)=SYS(JND1,KND2)+(REAL(IL)**2)*QQ(J0,K0)*RCATI(IG,JG,IL+1)/DELX
+                IF(IL.LE.NLF-3) THEN
+                  KND2=((IL+1)/2)*L4+(JG-1)*IELEM+K0
+                  SYS(JND1,KND2)=SYS(JND1,KND2)+REAL(IL*(IL+1))*QQ(J0,K0)*RCATI(IG,JG,IL+1)/DELX
+                ENDIF
+                IF(IL.GE.3) THEN
+                  KND2=((IL-3)/2)*L4+(JG-1)*IELEM+K0
+                  SYS(JND1,KND2)=SYS(JND1,KND2)+REAL((IL-1)*(IL-2))*QQ(J0,K0)*RCATI(IG,JG,IL-1)/DELX
+                  KND2=((IL-1)/2)*L4+(JG-1)*IELEM+K0
+                  SYS(JND1,KND2)=SYS(JND1,KND2)+(REAL(IL-1)**2)*QQ(J0,K0)*RCATI(IG,JG,IL-1)/DELX
+                ENDIF
+              ENDDO
+            ENDDO
+          ENDDO
+        ENDIF
+        DO J0=1,IELEM
+          JND1=((IL-1)/2)*L4+(IG-1)*IELEM+J0
+          GJXM=JXM(IG,1)
+          GJXP=JXP(IG,1)
+          IF(IL.EQ.3) THEN
+            GJXM=2.0*JXM(IG,1)+3.0*JXM(IG,2)
+            GJXP=2.0*JXP(IG,1)+3.0*JXP(IG,2)
+          ELSE IF(IL.EQ.5) THEN
+            GJXM=4.0*JXM(IG,2)+5.0*JXM(IG,3)
+            GJXP=4.0*JXP(IG,2)+5.0*JXP(IG,3)
+          ELSE IF(IL.EQ.7) THEN
+            GJXM=6.0*JXM(IG,3)+7.0*JXM(IG,4)
+            GJXP=6.0*JXP(IG,3)+7.0*JXP(IG,4)
+          ELSE IF(IL.GE.9) THEN
+            CALL XABORT('BRERTS: SPN ORDER NOT IMPLEMENTED(1).')
+          ENDIF
+          SYS(JND1,INX)=SYS(JND1,INX)-(V(1,J0)*GJXM+V(IELEM+1,J0)*GJXP)
+        ENDDO
+      ENDDO
+    ENDIF
+  ENDDO
+  !----
+  !  SOLVE A ONE-NODE PROBLEM.
+  !----
+  CALL ALSBD(L4*NLF/2,1,SYS,IER,L4*NLF/2)
+  IF(IER.NE.0) CALL XABORT('BRERTS: SINGULAR MATRIX(2).')
+  ALLOCATE(FUNKNO(IELEM,NGRP,NLF/2))
+  DO IL=1,NLF/2
+    DO IG=1,NGRP
+      DO J0=1,IELEM
+        FUNKNO(J0,IG,IL)=REAL(SYS((IL-1)*L4+(IG-1)*IELEM+J0,INX))
+      ENDDO
+    ENDDO
+  ENDDO
+  DEALLOCATE(SYS,RCATI,R,V)
+  !----
+  !  COMPUTE SURFACE FLUX GRADIENTS USING ODD PARITY EQUATIONS
+  !----
+  ALLOCATE(DXM(NGRP,NLF/2),DXP(NGRP,NLF/2))
+  IF(NLF.EQ.2) THEN
+    DXM(:NGRP,1)=MATMUL(RCAT(:,:,2),JXM(:,1))*DELX
+    DXP(:NGRP,1)=MATMUL(RCAT(:,:,2),JXP(:,1))*DELX
+  ELSE IF(NLF.EQ.4) THEN
+    DXM(:NGRP,2)=MATMUL(RCAT(:,:,4),JXM(:,2))*DELX/3.0
+    DXP(:NGRP,2)=MATMUL(RCAT(:,:,4),JXP(:,2))*DELX/3.0
+    DXM(:NGRP,1)=MATMUL(RCAT(:,:,2),JXM(:,1))*DELX-2.0*DXM(:NGRP,2)
+    DXP(:NGRP,1)=MATMUL(RCAT(:,:,2),JXP(:,1))*DELX-2.0*DXP(:NGRP,2)
+  ELSE IF(NLF.EQ.6) THEN
+    DXM(:NGRP,3)=MATMUL(RCAT(:,:,6),JXM(:,3))*DELX/5.0
+    DXP(:NGRP,3)=MATMUL(RCAT(:,:,6),JXP(:,3))*DELX/5.0
+    DXM(:NGRP,2)=MATMUL(RCAT(:,:,4),JXM(:,2))*DELX/3.0-4.0*DXM(:NGRP,3)/3.0
+    DXP(:NGRP,2)=MATMUL(RCAT(:,:,4),JXP(:,2))*DELX/3.0-4.0*DXP(:NGRP,3)/3.0
+    DXM(:NGRP,1)=MATMUL(RCAT(:,:,2),JXM(:,1))*DELX-2.0*DXM(:NGRP,2)
+    DXP(:NGRP,1)=MATMUL(RCAT(:,:,2),JXP(:,1))*DELX-2.0*DXP(:NGRP,2)
+  ELSE
+    CALL XABORT('BRERTS: SPN ORDER NOT IMPLEMENTED(2).')
+  ENDIF
+  !----
+  !  COMPUTE NODAL SURFACE FLUXES
+  !----
+  DENOM=REAL(IELEM*(IELEM+1))
+  DO IG=1,NGRP
+    DO IL=1,NLF/2
+      FHOMM_IG=0.0D0
+      FHOMP_IG=0.0D0
+      IF(ICOL.EQ.1) THEN
+        ! NEM relations
+        IF(IELEM.EQ.1) THEN
+          FHOMM_IG=FUNKNO(1,IG,IL)+((1./3.)*DXM(IG,IL)+(1./6.)*DXP(IG,IL))
+          FHOMP_IG=FUNKNO(1,IG,IL)-((1./6.)*DXM(IG,IL)+(1./3.)*DXP(IG,IL))
+        ELSE IF(IELEM.EQ.2) THEN
+          FHOMM_IG=FUNKNO(1,IG,IL)-(5.0*SQRT(3.)/6.)*FUNKNO(2,IG,IL)+((1./8.)*DXM(IG,IL)- &
+          & (1./24.)*DXP(IG,IL))
+          FHOMP_IG=FUNKNO(1,IG,IL)+(5.0*SQRT(3.)/6.)*FUNKNO(2,IG,IL)-(-(1./24.)*DXM(IG,IL)+ &
+          & (1./8.)*DXP(IG,IL))
+        ELSE IF(IELEM.EQ.3) THEN
+          FHOMM_IG=FUNKNO(1,IG,IL)-(5.0*SQRT(3.)/6.)*FUNKNO(2,IG,IL)+(7.0*SQRT(5.)/10.)* &
+          & FUNKNO(3,IG,IL)+((1./15.)*DXM(IG,IL)+(1./60.)*DXP(IG,IL))
+          FHOMP_IG=FUNKNO(1,IG,IL)+(5.0*SQRT(3.)/6.)*FUNKNO(2,IG,IL)+(7.0*SQRT(5.)/10.)* &
+          & FUNKNO(3,IG,IL)-((1./60.)*DXM(IG,IL)+(1./15.)*DXP(IG,IL))
+        ELSE
+          CALL XABORT('BRERTS: IELEM OVERFLOW.')
+        ENDIF
+      ELSE IF(ICOL.EQ.2) THEN
+        ! MCFD relations
+        FHOMM_IG=DXM(IG,IL)/DENOM
+        FHOMP_IG=-DXP(IG,IL)/DENOM
+        DO J0=1,IELEM
+          FP=SQRT(REAL(2*J0-1))*(1.0-REAL(J0*(J0-1))/DENOM)
+          FM=FP*(-1.0)**(J0-1)
+          FHOMM_IG=FHOMM_IG+FP*(-1.0)**(J0-1)*FUNKNO(J0,IG,IL)
+          FHOMP_IG=FHOMP_IG+FP*FUNKNO(J0,IG,IL)
+        ENDDO
+      ELSE IF(ICOL.EQ.3) THEN
+        IF(IELEM.EQ.1) THEN
+          FHOMM_IG=FUNKNO(1,IG,IL)+((1./4.)*DXM(IG,IL)+(1./4.)*DXP(IG,IL))
+          FHOMP_IG=FUNKNO(1,IG,IL)-((1./4.)*DXM(IG,IL)+(1./4.)*DXP(IG,IL))
+        ELSE IF(IELEM.EQ.2) THEN
+          FHOMM_IG=FUNKNO(1,IG,IL)-SQRT(3.0)*FUNKNO(2,IG,IL)+((1./12.)*DXM(IG,IL)- &
+          & (1./12.)*DXP(IG,IL))
+          FHOMP_IG=FUNKNO(1,IG,IL)+SQRT(3.0)*FUNKNO(2,IG,IL)-(-(1./12.)*DXM(IG,IL)+ &
+          & (1./12.)*DXP(IG,IL))
+        ELSE IF(IELEM.EQ.3) THEN
+          FHOMM_IG=FUNKNO(1,IG,IL)-(5.0*SQRT(3.)/6.0)*FUNKNO(2,IG,IL)+SQRT(5.)*FUNKNO(3,IG,IL)+ &
+          & ((1./24.)*DXM(IG,IL)+(1./24.)*DXP(IG,IL))
+          FHOMP_IG=FUNKNO(1,IG,IL)+(5.0*SQRT(3.)/6.0)*FUNKNO(2,IG,IL)+SQRT(5.)*FUNKNO(3,IG,IL)- &
+          & ((1./24.)*DXM(IG,IL)+(1./24.)*DXP(IG,IL))
+        ELSE
+          CALL XABORT('BRERTS: IELEM OVERFLOW.')
+        ENDIF
+      ELSE
+        CALL XABORT('BRERTS: ICOL OVERFLOW.')
+      ENDIF
+      FHOMM(IG,IL)=REAL(FHOMM_IG)
+      FHOMP(IG,IL)=REAL(FHOMP_IG)
+    ENDDO
+    IF(IMPX.GT.0) THEN
+      DO IL=1,NLF/2
+        WRITE(6,'(14H BRERTS: ORDER,I2,38H RAVIART-THOMAS FLUX UNKNOWNS IN GROUP,I4,1H:/ &
+        & (1P,12E12.4))') 2*IL-2,IG,FUNKNO(:IELEM,IG,IL)
+        WRITE(6,'(14H BRERTS: ORDER,I2,40H RAVIART-THOMAS BOUNDARY FLUXES IN GROUP,I4,1H:/ &
+        & (1P,12E12.4))') 2*IL-2,IG,FHOMM(IG,IL),FHOMP(IG,IL)
+      ENDDO
+    ENDIF
+  ENDDO
+  DEALLOCATE(DXP,DXM,FUNKNO)
+END SUBROUTINE BRERTS