From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/SALMUS.f90 | 282 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 282 insertions(+) create mode 100644 Dragon/src/SALMUS.f90 (limited to 'Dragon/src/SALMUS.f90') 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 -- cgit v1.2.3