Initial commit from Polytechnique Montreal

1 files changed, 282 insertions, 0 deletions

diff --git a/Dragon/src/SALMUS.f90 b/Dragon/src/SALMUS.f90
new file mode 100644
index 0000000..df30a61
--- /dev/null
+++ b/Dragon/src/SALMUS.f90
@@ -0,0 +1,282 @@
+!---------------------------------------------------------------------
+!
+!Purpose:
+! To analyze and track a geometry data structure using the Sanchez
+! algorithm for the multicell surfacic approximation.
+!
+!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
+! FGEO    unit file number of the surfacic file in read only mode.
+! IPTRK   pointer to the TRACKING data structure in creation mode.
+! IFTRK   pointer to the TRACKING file in creation mode.
+! RCUTOF  minimum distance between two surfacic elements.
+! IPRINT  print level.
+! NBSLIN  maximum number of segments in a single tracking line.
+!         computed by default in SALTCG but limited to 100000
+!         elements. This default value can be bypassed using
+!         keyword NBSLIN.
+!
+!Parameters: output
+! GG      geometry basic information.
+!
+!-----------------------------------------------------------------------
+!
+SUBROUTINE SALMUS(FGEO ,IPTRK, IFTRK, RCUTOF, IPRINT, NBSLIN, GG)
+  USE GANLIB
+  USE PRECISION_AND_KINDS, ONLY : PDB
+  USE SAL_GEOMETRY_TYPES,  ONLY : T_G_BASIC,LMERGM,TYPGEO,NBFOLD,F_GEO
+  USE SAL_TRACKING_TYPES,  ONLY : PRTIND,EPS1
+  USE SAL_GEOMETRY_MOD,    ONLY : SAL100
+  IMPLICIT     NONE
+  !----
+  !  Subroutine arguments
+  !----
+  TYPE(C_PTR) IPTRK
+  INTEGER  FGEO,IFTRK,IPRINT,NBSLIN
+  REAL(PDB) RCUTOF
+  TYPE(T_G_BASIC) :: GG
+  !----
+  !  Local parameters
+  !----
+  INTEGER      IOUT
+  PARAMETER   (IOUT=6)
+  INTEGER      NSTATE
+  PARAMETER   (NSTATE=40)
+  !----
+  !  Local variables
+  !----
+  TYPE(C_PTR)  JPTRK,KPTRK
+  INTEGER      ISTATT(NSTATE),I,J,I1,I2,J1,J2,IMACRO,J3,NMIX,IMIX,NREG,OK, &
+               & NINST,ILONG,ITYLCM
+  REAL         RSTATT(NSTATE)
+  LOGICAL      LGINF
+  CHARACTER(LEN=131) HSMG
+  !----
+  !  Allocatable arrays
+  !----
+  INTEGER, ALLOCATABLE, DIMENSION(:) :: NSURF_MACRO,NBNODE_MACRO,ITAB, &
+  & NMC_NODE,NMC_SURF,IFR,MIX,ICNT,MERGE_MACRO,SURF_MACRO,PERIM,IDMA,IMAC
+  REAL, ALLOCATABLE, DIMENSION(:) :: VOLUME,ALB,DVX,SUR,GALBED
+  REAL(PDB), ALLOCATABLE, DIMENSION(:) :: SURF2
+  !----
+  !  Read the surfacic file and fill GG object
+  !----
+  PRTIND=IPRINT
+  F_GEO=FGEO
+  EPS1=1.E-5_PDB
+  IF(RCUTOF>0._PDB) THEN
+    EPS1=RCUTOF
+    IF(PRTIND>0) WRITE(*,*) "SALMUS: set eps1 to ",EPS1
+  ENDIF
+  CALL SAL100(GG)
+  !------------
+  IF(IPRINT > 0) WRITE(IOUT,'(/" SALMUS: TYPGEO=",I5," NBFOLD=",I5," NMACRO=",I5)') &
+  & TYPGEO,NBFOLD,GG%NB_MACRO
+  !----
+  ! Perform optional MERGE MIX
+  !----
+  IF(LMERGM) THEN
+    GG%NUM_MERGE(:)=GG%MED(:)
+    DO I=1,MAXVAL(GG%NUM_MERGE(:))
+      1000 IF(I.GT.MAXVAL(GG%NUM_MERGE(:))) EXIT
+      DO J=1,GG%NB_NODE
+        IF(GG%NUM_MERGE(J).EQ.I) GO TO 2000
+      ENDDO
+      DO J=1,GG%NB_NODE
+        IF(GG%NUM_MERGE(J).GE.I) GG%NUM_MERGE(J)=GG%NUM_MERGE(J)-1
+      ENDDO
+      GO TO 1000
+      2000 CONTINUE
+    ENDDO
+  ENDIF
+  !----
+  ! store the STATE VECTOR for the global geometry
+  !----
+  CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATT)
+  CALL LCMGET(IPTRK,'EXCELTRACKOP',RSTATT)
+  NREG=MAXVAL(GG%NUM_MERGE)
+  ISTATT(1) = NREG ! number of regions
+  ISTATT(2) = NREG ! number of flux unknowns
+  ISTATT(7) = 5 ! set the multicell surfacic approximation
+  LGINF = .TRUE.
+  DO I=1,GG%NBBCDA
+    LGINF = LGINF .AND. (GG%BCDATAREAD(I)%BCDATA(6) == 1._PDB)
+  ENDDO
+  ISTATT(3)=1
+  IF(.NOT.LGINF) ISTATT(3)=0 ! reset the leakage flag
+  ISTATT(4) = MAXVAL(GG%MED(1:GG%NB_NODE)) ! maximum number of mixture
+  CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,ISTATT)
+  !
+  ! fill-in medium number per region
+  ALLOCATE(ITAB(NREG),VOLUME(NREG), STAT =OK)
+  IF(OK /= 0) CALL XABORT('SALMUS: failure to allocate integer ITAB')
+  ! fill in MATCOD
+  DO J=1,GG%NB_NODE
+    ITAB(GG%NUM_MERGE(J)) = GG%MED(J)
+  ENDDO
+  CALL LCMPUT(IPTRK,'MATCOD',NREG,1,ITAB(1:NREG) )
+  ! fill-in KEYFLX per region
+  ITAB(:NREG) = (/ (I,I=1,NREG) /)
+  CALL LCMPUT(IPTRK,'MERGE',NREG,1,ITAB)
+  CALL LCMPUT(IPTRK,'KEYFLX',NREG,1,ITAB)
+  ! fill-in volumes per region
+  VOLUME(:NREG) =0.
+  DO I=1,GG%NB_NODE
+    VOLUME(GG%NUM_MERGE(I)) = VOLUME(GG%NUM_MERGE(I)) + REAL(GG%VOL_NODE(I))
+  ENDDO
+  CALL LCMPUT(IPTRK,'VOLUME',NREG,2,VOLUME)
+  IF(IPRINT .GT. 5) THEN
+    I1=1
+    DO I=1,(NREG-1)/8+1
+      I2=I1+7
+      IF(I2.GT.NREG) I2=NREG
+      WRITE (IOUT,20) (J,J=I1,I2)
+      DO J=1,GG%NB_NODE
+        ITAB(GG%NUM_MERGE(J)) = GG%MED(J)
+      ENDDO
+      WRITE (IOUT,30) (ITAB(J),J=I1,I2)
+      WRITE (IOUT,40) (VOLUME(J),J=I1,I2)
+      I1=I1+8
+    ENDDO
+  ENDIF
+  DEALLOCATE(VOLUME,ITAB)
+  !----
+  ! Extract the surfacic elements belonging to each macro geometry and
+  ! perform tracking
+  !----
+  ALLOCATE(NSURF_MACRO(GG%NB_MACRO),NBNODE_MACRO(GG%NB_MACRO))
+  ALLOCATE(NMC_NODE(GG%NB_MACRO+1),NMC_SURF(GG%NB_MACRO+1))
+  NSURF_MACRO(:GG%NB_MACRO) = 0
+  NBNODE_MACRO(:GG%NB_MACRO) = 0
+  JPTRK=LCMLID(IPTRK,'MACRO-TRACK',GG%NB_MACRO)
+  NMC_NODE(1)=0
+  NMC_SURF(1)=0
+  DO IMACRO=1,GG%NB_MACRO
+    KPTRK=LCMDIL(JPTRK,IMACRO)
+    CALL LCMPUT(KPTRK,'STATE-VECTOR',NSTATE,1,ISTATT)
+    CALL LCMPUT(KPTRK,'EXCELTRACKOP',NSTATE,2,RSTATT)
+    CALL MUSACG(KPTRK,IFTRK,IPRINT,IMACRO,NBSLIN,RCUTOF,GG,LGINF, &
+    & NBNODE_MACRO(IMACRO),NSURF_MACRO(IMACRO))
+    NMC_NODE(IMACRO+1)=NMC_NODE(IMACRO)+NBNODE_MACRO(IMACRO)
+    NMC_SURF(IMACRO+1)=NMC_SURF(IMACRO)+NSURF_MACRO(IMACRO)
+  ENDDO
+  CALL LCMPUT(IPTRK,'NMC_NODE',GG%NB_MACRO+1,1,NMC_NODE)
+  CALL LCMPUT(IPTRK,'NMC_SURF',GG%NB_MACRO+1,1,NMC_SURF)
+  !----
+  !  Create connectivity data
+  !----
+  NMIX=NMC_SURF(GG%NB_MACRO+1)
+  ALLOCATE(IFR(NMIX),ALB(NMIX),MIX(NMIX),DVX(NMIX),SUR(NMIX),ICNT(NMIX),IDMA(NMIX))
+  ALLOCATE(IMAC(NREG))
+  J1=0
+  NMIX=0
+  IMIX=0
+  DO IMACRO=1,GG%NB_MACRO
+    KPTRK=LCMDIL(JPTRK,IMACRO)
+    J2=NBNODE_MACRO(IMACRO)
+    J3=NSURF_MACRO(IMACRO)
+    ALLOCATE(MERGE_MACRO(J2),SURF_MACRO(J3),SURF2(J3))
+    CALL LCMSIX(KPTRK,'SURFACIC_TMP',1)
+      CALL LCMGET(KPTRK,'MERGE_MACRO',MERGE_MACRO)
+      CALL LCMGET(KPTRK,'SURF_MACRO',SURF_MACRO)
+    CALL LCMSIX(KPTRK,' ',2)
+    CALL LCMSIX(KPTRK,'GEOMETRY',1)
+      CALL LCMLEN(KPTRK,'SURF2',ILONG,ITYLCM)
+      IF(ILONG.NE.J3) CALL XABORT('SALMUS: wrong number of surfaces')
+      CALL LCMGET(KPTRK,'SURF2',SURF2)
+    CALL LCMSIX(KPTRK,' ',2)
+    IF(J1+J2.GT.NREG) CALL XABORT('SALMUS: NREG overflow')
+    IMAC(J1+1:J1+J2)=MERGE_MACRO(:J2)
+    ICNT(NMIX+1:NMIX+J3)=SURF_MACRO(:J3)
+    DEALLOCATE(SURF_MACRO,MERGE_MACRO)
+    CALL LCMLEN(KPTRK,'ALBEDO',ILONG,ITYLCM)
+    ALLOCATE(GALBED(ILONG),PERIM(J3))
+    CALL LCMGET(KPTRK,'ALBEDO',GALBED)
+    CALL LCMGET(KPTRK,'PERIM_SURF',PERIM)
+    DO I=1,J3
+      IF(PERIM(I).GT.ILONG) THEN
+        WRITE(HSMG,'(51H SALMUS: inconsistent albedo info in macro geometry,I5)') IMACRO
+        CALL XABORT(HSMG)
+      ENDIF
+      ALB(NMIX+I)=GALBED(PERIM(I))
+      SUR(NMIX+I)=REAL(SURF2(I))
+    ENDDO
+    DEALLOCATE(PERIM,GALBED,SURF2)
+    OUT1: DO I=NMIX+1,NMIX+J3
+      DO J=NMIX+1,I-1
+        IF(ICNT(I).EQ.ICNT(J)) THEN
+          MIX(I)=MIX(J)
+          CYCLE OUT1
+        ENDIF
+      ENDDO
+      IMIX=IMIX+1
+      MIX(I)=IMIX
+    ENDDO OUT1
+    DO I=NMIX+1,NMIX+J3
+      NINST=COUNT(MIX(NMIX+1:NMIX+J3) == MIX(I))
+      DVX(I)=1.0/REAL(NINST)
+      IDMA(I)=IMACRO
+    ENDDO
+    J1=J1+J2
+    NMIX=NMIX+J3
+  ENDDO
+  CALL LCMPUT(IPTRK,'MERGE_MACRO',NREG,1,IMAC)
+  DEALLOCATE(IMAC)
+  OUT2: DO I=1,NMIX
+    DO J=1,NMIX
+      IF(IDMA(J).EQ.IDMA(I)) CYCLE
+      IF(ICNT(J).EQ.ICNT(I)) THEN
+        IFR(I)=MIX(J)
+        CYCLE OUT2
+      ENDIF
+    ENDDO
+    IFR(I)=MIX(I)
+  ENDDO OUT2
+  DEALLOCATE(IDMA)
+  ISTATT(24)=GG%NB_MACRO
+  ISTATT(28)=MAXVAL(IFR) ! number of current unknowns
+  ISTATT(29)=NMIX ! number of perimeter elements
+  CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,ISTATT)
+  CALL LCMPUT(IPTRK,'IFR',NMIX,1,IFR)
+  CALL LCMPUT(IPTRK,'ALB',NMIX,2,ALB)
+  CALL LCMPUT(IPTRK,'MIX',NMIX,1,MIX)
+  CALL LCMPUT(IPTRK,'DVX',NMIX,2,DVX)
+  CALL LCMPUT(IPTRK,'SUR',NMIX,2,SUR)
+  IF(IPRINT .GT. 1) THEN
+    NMIX=0
+    DO IMACRO=1,GG%NB_MACRO
+      J3=NSURF_MACRO(IMACRO)
+      WRITE (IOUT,50) IMACRO,(ICNT(I),I=NMIX+1,NMIX+J3)
+      WRITE (IOUT,60) (MIX(I),IFR(I),I=NMIX+1,NMIX+J3)
+      WRITE (IOUT,70) (ALB(I),I=NMIX+1,NMIX+J3)
+      WRITE (IOUT,80) (DVX(I),I=NMIX+1,NMIX+J3)
+      WRITE (IOUT,90) (SUR(I),I=NMIX+1,NMIX+J3)
+      WRITE (IOUT,100) ('----------------',I=1,MIN(8,J3))
+      NMIX=NMIX+J3
+    ENDDO
+  ENDIF
+  DEALLOCATE(ICNT,SUR,DVX,MIX,ALB,IFR)
+  DEALLOCATE(NMC_SURF,NMC_NODE)
+  DEALLOCATE(NBNODE_MACRO,NSURF_MACRO)
+  !----
+  !  Formats
+  !----
+   20 FORMAT (/9H  REGION:,8(I8,6X,1HI))
+   30 FORMAT (9H MIXTURE:,8(I8,6X,1HI))
+   40 FORMAT (9H  VOLUME:,1P,8(E13.6,2H I))
+   50 FORMAT (6H MACRO,I6.6/9H ELEMENT:,8(3H  S,I6.6,6X,1HI,:)/(9X,8(3H  S,I6.6,6X,1HI,:)))
+   60 FORMAT (9H  IN/OUT:,8(I6,2H /,I5,3H  I,:)/(9X,8(I6,2H /,I5,3H  I,:)))
+   70 FORMAT (9H  ALBEDO:,1P,8(E13.5,3H  I,:)/(9X,8(E13.5,3H  I,:)))
+   80 FORMAT (9H     DVX:,1P,8(E13.5,3H  I,:)/(9X,8(E13.5,3H  I,:)))
+   90 FORMAT (9H     SUR:,1P,8(E13.5,3H  I,:)/(9X,8(E13.5,3H  I,:)))
+  100 FORMAT (9H --------,8(A16))
+END SUBROUTINE SALMUS