Initial commit from Polytechnique Montreal

1 files changed, 180 insertions, 0 deletions

diff --git a/Dragon/src/BRERTD.f90 b/Dragon/src/BRERTD.f90
new file mode 100644
index 0000000..9a62fca
--- /dev/null
+++ b/Dragon/src/BRERTD.f90
@@ -0,0 +1,180 @@
+SUBROUTINE BRERTD(IELEM,ICOL,NGRP,DELX,DIFF,SIGT,SUNXS,JXM,JXP,IMPX, &
+& FHOMM,FHOMP)
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! Compute the Raviart-Thomas boundary fluxes for a single node.
+!
+!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.
+! DELX    node width along X-axis.
+! DIFF    diffusion coefficient array (cm).
+! SIGT    total XS (cm^-1).
+! SUNXS   production (fission + scattering) cross sections. The second
+!         dimension is for primary neutrons.
+! 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,IMPX
+  REAL, INTENT(IN) :: DELX
+  REAL, DIMENSION(NGRP), INTENT(IN) :: DIFF,JXM,JXP
+  REAL, DIMENSION(NGRP), INTENT(IN) :: SIGT
+  REAL, DIMENSION(NGRP,NGRP), INTENT(IN) :: SUNXS
+  REAL, DIMENSION(NGRP), INTENT(OUT) :: FHOMM,FHOMP
+  !----
+  !  LOCAL VARIABLES
+  !----
+  TYPE(C_PTR) IPTRK
+  REAL QQ(5,5)
+  !----
+  !  ALLOCATABLE ARRAYS
+  !----
+  REAL, DIMENSION(:,:), ALLOCATABLE :: R,V,SYS,FUNKNO
+  !----
+  !  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)
+  !----
+  !  LEAKAGE-REMOVAL SYSTEM MATRIX ASSEMBLY.
+  !----
+  ALLOCATE(SYS(IELEM*NGRP,IELEM*NGRP+1))
+  SYS(:IELEM*NGRP,:IELEM*NGRP+1)=0.0
+  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
+  INX=IELEM*NGRP+1
+  DO IG=1,NGRP
+    DO J0=1,IELEM
+      JND1=(IG-1)*IELEM+J0
+      SYS(JND1,JND1)=SYS(JND1,JND1)+DELX*SIGT(IG)
+      DO K0=1,IELEM
+        IF(QQ(J0,K0).EQ.0.0) CYCLE
+        KND1=(IG-1)*IELEM+K0
+        SYS(JND1,KND1)=SYS(JND1,KND1)+QQ(J0,K0)*DIFF(IG)/DELX
+      ENDDO
+      SYS(JND1,INX)=SYS(JND1,INX)-V(1,J0)*JXM(IG)-V(IELEM+1,J0)*JXP(IG)
+    ENDDO
+  ENDDO
+  !----
+  !  SOURCE SYSTEM MATRIX ASSEMBLY.
+  !----
+  DO IG=1,NGRP
+    DO JG=1,NGRP
+      DO J0=1,IELEM
+        JND1=(IG-1)*IELEM+J0
+        JND2=(JG-1)*IELEM+J0
+        SYS(JND1,JND2)=SYS(JND1,JND2)-DELX*SUNXS(IG,JG) ! IG <-- JG
+      ENDDO
+    ENDDO
+  ENDDO
+  !----
+  !  SOLVE A ONE-NODE PROBLEM.
+  !----
+  CALL ALSB(IELEM*NGRP,1,SYS,IER,IELEM*NGRP)
+  IF(IER.NE.0) CALL XABORT('BRERTD: SINGULAR MATRIX.')
+  ALLOCATE(FUNKNO(IELEM,NGRP))
+  DO IG=1,NGRP
+    DO J0=1,IELEM
+      FUNKNO(J0,IG)=SYS((IG-1)*IELEM+J0,IELEM*NGRP+1)
+    ENDDO
+  ENDDO
+  DEALLOCATE(SYS,R,V)
+  !----
+  !  COMPUTE NODAL SURFACE FLUXES
+  !----
+  DENOM=REAL(IELEM*(IELEM+1))
+  FHOMM(:NGRP)=0.0
+  FHOMP(:NGRP)=0.0
+  DO IG=1,NGRP
+    IF(ICOL.EQ.1) THEN
+      ! NEM relations
+      IF(IELEM.EQ.1) THEN
+        FHOMM(IG)=FUNKNO(1,IG)+((1./3.)*JXM(IG)+(1./6.)*JXP(IG))*DELX/DIFF(IG)
+        FHOMP(IG)=FUNKNO(1,IG)-((1./6.)*JXM(IG)+(1./3.)*JXP(IG))*DELX/DIFF(IG)
+      ELSE IF(IELEM.EQ.2) THEN
+        FHOMM(IG)=FUNKNO(1,IG)-(5.0*SQRT(3.0)/6.0)*FUNKNO(2,IG)+((1./8.)*JXM(IG)- &
+        & (1./24.)*JXP(IG))*DELX/DIFF(IG)
+        FHOMP(IG)=FUNKNO(1,IG)+(5.0*SQRT(3.0)/6.0)*FUNKNO(2,IG)-(-(1./24.)*JXM(IG)+ &
+        & (1./8.)*JXP(IG))*DELX/DIFF(IG)
+      ELSE IF(IELEM.EQ.3) THEN
+        FHOMM(IG)=FUNKNO(1,IG)-(5.0*SQRT(3.0)/6.0)*FUNKNO(2,IG)+(7.0*SQRT(5.0)/10.0)*FUNKNO(3,IG)+ &
+        & ((1./15.)*JXM(IG)+(1./60.)*JXP(IG))*DELX/DIFF(IG)
+        FHOMP(IG)=FUNKNO(1,IG)+(5.0*SQRT(3.0)/6.0)*FUNKNO(2,IG)+(7.0*SQRT(5.0)/10.0)*FUNKNO(3,IG)- &
+        & ((1./60.)*JXM(IG)+(1./15.)*JXP(IG))*DELX/DIFF(IG)
+      ELSE
+        CALL XABORT('BRERTD: IELEM OVERFLOW.')
+      ENDIF
+    ELSE IF(ICOL.EQ.2) THEN
+      ! MCFD relations
+      FHOMM(IG)=JXM(IG)*DELX/(DIFF(IG)*DENOM)
+      FHOMP(IG)=-JXP(IG)*DELX/(DIFF(IG)*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)
+        FHOMP(IG)=FHOMP(IG)+FP*FUNKNO(J0,IG)
+      ENDDO
+    ELSE IF(ICOL.EQ.3) THEN
+      IF(IELEM.EQ.1) THEN
+        FHOMM(IG)=FUNKNO(1,IG)+((1./4.)*JXM(IG)+(1./4.)*JXP(IG))*DELX/DIFF(IG)
+        FHOMP(IG)=FUNKNO(1,IG)-((1./4.)*JXM(IG)+(1./4.)*JXP(IG))*DELX/DIFF(IG)
+      ELSE IF(IELEM.EQ.2) THEN
+        FHOMM(IG)=FUNKNO(1,IG)-SQRT(3.0)*FUNKNO(2,IG)+((1./12.)*JXM(IG)- &
+        & (1./12.)*JXP(IG))*DELX/DIFF(IG)
+        FHOMP(IG)=FUNKNO(1,IG)+SQRT(3.0)*FUNKNO(2,IG)-(-(1./12.)*JXM(IG)+ &
+        & (1./12.)*JXP(IG))*DELX/DIFF(IG)
+      ELSE IF(IELEM.EQ.3) THEN
+        FHOMM(IG)=FUNKNO(1,IG)-(5.0*SQRT(3.0)/6.0)*FUNKNO(2,IG)+SQRT(5.0)*FUNKNO(3,IG)+ &
+        & ((1./24.)*JXM(IG)+(1./24.)*JXP(IG))*DELX/DIFF(IG)
+        FHOMP(IG)=FUNKNO(1,IG)+(5.0*SQRT(3.0)/6.0)*FUNKNO(2,IG)+SQRT(5.0)*FUNKNO(3,IG)- &
+        & ((1./24.)*JXM(IG)+(1./24.)*JXP(IG))*DELX/DIFF(IG)
+      ELSE
+        CALL XABORT('BRERTD: IELEM OVERFLOW.')
+      ENDIF
+    ELSE
+      CALL XABORT('BRERTD: ICOL OVERFLOW.')
+    ENDIF
+    IF(IMPX.GT.0) THEN
+      WRITE(6,'(46H BRERTD: RAVIART-THOMAS FLUX UNKNOWNS IN GROUP,I4,1H:/(1P,12E12.4))') &
+      & IG,FUNKNO(:IELEM,IG)
+      WRITE(6,'(48H BRERTD: RAVIART-THOMAS BOUNDARY FLUXES IN GROUP,I4,1H:/(1P,12E12.4))') &
+      & IG,FHOMM(IG),FHOMP(IG)
+    ENDIF
+  ENDDO
+  DEALLOCATE(FUNKNO)
+END SUBROUTINE BRERTD