Initial commit from Polytechnique Montreal

1 files changed, 487 insertions, 0 deletions

diff --git a/Donjon/src/NAPGEO.f b/Donjon/src/NAPGEO.f
new file mode 100644
index 0000000..da231e6
--- /dev/null
+++ b/Donjon/src/NAPGEO.f
@@ -0,0 +1,487 @@
+*DECK NAPGEO

+      SUBROUTINE NAPGEO(IPGNW,IPGOD,IPCPO,NSTATE)

+*

+*-----------------------------------------------------------------------

+*

+*Purpose:

+* Split geometry from homogeneous to heterogeneous assemblies

+*

+*Copyright:

+* Copyright (C) 2014 Ecole Polytechnique de Montreal.

+*

+*Author(s): 

+* R. Chambon

+*

+*Parameters: input/output

+* IPGNW   LCM object address of heterogeneous assembly Geometry.

+* IPGOD   LCM object address of homogeneous assembly Geometry.

+* IPCPO   LCM object address of Multicompo.

+* NSTATE  length of the state vector

+*

+*-----------------------------------------------------------------------

+*

+      USE GANLIB

+      IMPLICIT NONE

+*----

+*  SUBROUTINE ARGUMENTS

+*----

+      INTEGER      NSTATE

+      TYPE(C_PTR)  IPGNW,IPGOD,IPCPO

+*----

+*  LOCAL VARIABLES

+*----

+      INTEGER IOUT,MAXLIN

+      REAL REPS

+      PARAMETER (REPS=1.0E-4,IOUT=6,MAXLIN=50)

+      TYPE(C_PTR) JPGEO,KPGEO,JPCPO,KPCPO

+      INTEGER INDIC,NITMA,LENGTH

+      CHARACTER TEXT*12

+      REAL FLOT

+      DOUBLE PRECISION DFLOT

+      INTEGER ISTATE(NSTATE),IMPX,NCODE(6),ICODE(6),ITYPGP,STYPP,

+     1  KCHAR(3)

+      REAL ZCODE(6)

+      INTEGER NXP,NYP,NREGP,NMIXP,NCOMLI,NXD,NYD,NZD,NREGD,NMIXD,NXF,

+     1  NYF,NZF,NREGF,NMIXF,NXA,NYA,NMIXA,NXPTMP,NYPTMP,NMIXD2,NASS

+      CHARACTER DIRHET*12,COMMEN(MAXLIN)*80,HMSG*131

+      INTEGER I,J,K,L,IP,JP,JF,IFBEG,JFBEG,LMIX,IASS,IZ,IZT,JM,JN,IM

+      LOGICAL LSPX,LSPY,LSPZ,LPOS,LMGEO

+      INTEGER, ALLOCATABLE, DIMENSION(:) :: MIXA,NBAX,AZONE,IBAX

+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: MIXP

+      INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: MIXD,MIXF

+      REAL, ALLOCATABLE, DIMENSION(:) :: MXP,MYP,MXD,MYD,MZD,MXF,

+     1  MYF

+      INTEGER, ALLOCATABLE, DIMENSION(:) :: SXP,SYP,SXD,SYD,SZD,

+     1  SXF,SYF,AXD,AYD

+

+      IMPX=0

+      LSPX=.FALSE.

+      LSPY=.FALSE.

+      LSPZ=.FALSE.

+      LMGEO=.FALSE.

+C Read mandatory inputs

+      CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)

+      IF(INDIC.NE.3) CALL XABORT('NAPGEO: character data expected.')

+      IF(TEXT.EQ.'EDIT') THEN

+        CALL REDGET(INDIC,IMPX,FLOT,TEXT,DFLOT)

+        IF(INDIC.NE.1) CALL XABORT('NAPGEO: integer data expected.')

+        CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)

+        IF(INDIC.NE.3) CALL XABORT('NAPGEO: character data expected.')

+      ENDIF

+      IF(TEXT.NE.'DIRGEO') CALL XABORT('NAPGEO: ''DIRGEO'' '//

+     1  'expected.')

+      CALL REDGET(INDIC,NITMA,FLOT,DIRHET,DFLOT)

+      IF(INDIC.NE.3) CALL XABORT('NAPGEO: character data expected.')

+      CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)

+      IF(TEXT.EQ.'MACGEO') THEN

+        LMGEO=.TRUE.

+        CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)

+      ENDIF

+      IF(TEXT.NE.'MIXASS') THEN

+        CALL XABORT('NAPGEO: ''MIXASS'' expected.')

+      ENDIF

+      CALL REDGET(INDIC,NMIXA,FLOT,DIRHET,DFLOT)

+      IF(INDIC.NE.1)CALL XABORT('@NAPGEO: integer data expected.')

+      ALLOCATE(MIXA(NMIXA*2))

+      DO I=1,NMIXA

+        CALL REDGET(INDIC,MIXA(I),FLOT,DIRHET,DFLOT)

+        IF(INDIC.NE.1)CALL XABORT('@NAPGEO: integer data expected.')

+      ENDDO

+      CALL LCMSIX(IPCPO,DIRHET,1)

+

+*----

+*  RECOVER TABLE-OF-CONTENT INFORMATION FOR THE COMPO.

+*----

+      CALL LCMGET(IPCPO,'STATE-VECTOR',ISTATE)

+      NCOMLI=ISTATE(10)

+      CALL LCMGTC(IPCPO,'COMMENT',80,NCOMLI,COMMEN)

+      IF(IMPX.GT.10)WRITE(IOUT,'(1X,A)') (COMMEN(I),I=1,NCOMLI)

+* Get Geometry from calculation #1

+* for pin by pin geometries

+      IF(LMGEO) THEN

+        JPCPO=LCMGID(IPCPO,'MIXTURES')

+        KPCPO=LCMGIL(JPCPO,1)

+        JPGEO=LCMGID(KPCPO,'CALCULATIONS')

+        KPGEO=LCMGIL(JPGEO,1)

+        CALL LCMSIX(KPGEO,'MACROLIB    ',1)

+        CALL LCMSIX(KPGEO,'GFF         ',1)

+        CALL LCMSIX(KPGEO,'GFF-GEOM    ',1)

+      ELSE

+* for heterogeneous geometries

+        JPGEO=LCMGID(IPCPO,'GEOMETRIES')

+        KPGEO=LCMGIL(JPGEO,1)

+      ENDIF

+C get dimension in geometry from L_MULTICOMPO

+      if(impx.ge.100)write(6,*) 'debug: get ISTATE multicompo Geometry'

+      CALL LCMGET(KPGEO,'STATE-VECTOR',ISTATE)

+      ITYPGP=ISTATE(1)

+      STYPP=ISTATE(11)

+      IF(ITYPGP.NE.5) CALL XABORT('NAPGEO: CAR2D geometry type '

+     1 //'expected in L_MULTICOMPO.')

+      IF(STYPP.NE.0) CALL XABORT('NAPGEO: No split in geometry '

+     1 //'expected in L_MULTICOMPO.')

+      NXP=ISTATE(3)

+      NYP=ISTATE(4)

+      NREGP=ISTATE(6)

+      NMIXP=ISTATE(7)

+      ALLOCATE(MXP(NXP+1),MYP(NYP+1))

+      ALLOCATE(SXP(NXP),SYP(NYP))

+      ALLOCATE(MIXP(NXP,NYP))

+      CALL LCMGET(KPGEO,'MESHX',MXP)

+      CALL LCMGET(KPGEO,'MESHY',MYP)

+      CALL LCMGET(KPGEO,'MIX',MIXP)

+      SXP(:NXP)=1

+      SYP(:NYP)=1

+C get dimension in homogeneous assembly core Geometry

+      if(impx.ge.100)write(6,*) 'debug: ISTATE homog. ass. core Geom.'

+      CALL LCMGET(IPGOD,'SIGNATURE',KCHAR)

+      CALL LCMGET(IPGOD,'STATE-VECTOR',ISTATE)

+      ITYPGP=ISTATE(1)

+      NXD=ISTATE(3)

+      NYD=ISTATE(4)

+      NZD=ISTATE(5)

+      NREGD=ISTATE(6)

+      NMIXD=ISTATE(7)

+      NMIXD2=NMIXD

+      ALLOCATE(MXD(NXD+1),MYD(NYD+1))

+      ALLOCATE(MZD(NZD+1))

+      ALLOCATE(SXD(NXD),SYD(NYD))

+      ALLOCATE(SZD(NZD))

+      ALLOCATE(AXD(NXD),AYD(NYD))

+      ALLOCATE(MIXD(NXD,NYD,NZD))

+      CALL LCMGET(IPGOD,'MESHX',MXD)

+      CALL LCMGET(IPGOD,'MESHY',MYD)

+      CALL LCMGET(IPGOD,'MESHZ',MZD)

+      CALL LCMGET(IPGOD,'MIX',MIXD)

+      AXD(:NXD)=0

+      AYD(:NYD)=0

+      CALL LCMLEN(IPGOD,'SPLITX',LENGTH,INDIC)

+      IF(LENGTH.NE.0) THEN

+        CALL LCMGET(IPGOD,'SPLITX',SXD)

+        LSPX=.TRUE.

+      ELSE

+        SXD(:NXD)=1

+      ENDIF

+      CALL LCMLEN(IPGOD,'SPLITY',LENGTH,INDIC)

+      IF(LENGTH.NE.0) THEN

+        CALL LCMGET(IPGOD,'SPLITY',SYD)

+        LSPY=.TRUE.

+      ELSE

+        SYD(:NYD)=1

+      ENDIF

+      CALL LCMLEN(IPGOD,'SPLITZ',LENGTH,INDIC)

+      IF(LENGTH.NE.0) THEN

+        CALL LCMGET(IPGOD,'SPLITZ',SZD)

+        LSPZ=.TRUE.

+      ELSE

+        SZD(:NZD)=1

+      ENDIF

+      CALL LCMGET(IPGOD,'NCODE',NCODE)

+      CALL LCMGET(IPGOD,'ZCODE',ZCODE)

+      CALL LCMGET(IPGOD,'ICODE',ICODE)

+C get assembly mixture in homogeneous core geometry

+      if(impx.ge.100)write(6,*) 'debug: get assembly mixture'

+      DO 40 K=1,NZD

+      DO 30 J=1,NYD

+      DO 20 I=1,NXD

+      DO 10 L=1,NMIXA

+         IF(MIXA(L).EQ.MIXD(I,J,K)) THEN

+           AXD(I)=1

+           AYD(J)=1

+           GOTO 20

+         ENDIF

+   10 CONTINUE

+   20 CONTINUE

+   30 CONTINUE

+   40 CONTINUE

+      if(impx.ge.5) then

+        write(6,*) 'Original mesh corresponding to assemblies'

+        write(6,*) 'X direction: AXD(1 : NXD)=',(AXD(I),I=1,NXD)

+        write(6,*) 'Y direction: AYD(1 : NYD)=',(AYD(I),I=1,NYD)

+      endif

+C specify splitting in heterogeneous assembly geometry

+   50 CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)

+      IF(INDIC.NE.3) CALL XABORT('@NAPGEO: character data expected 1.')

+      IF(TEXT.EQ.'SPLITX-ASS') THEN

+        DO I=1,NXP

+          CALL REDGET(INDIC,SXP(I),FLOT,TEXT,DFLOT)

+          IF(INDIC.NE.1) THEN

+            WRITE(HMSG,*) '@NAPGEO: integer number expected'

+     1  //' for SPLITX-ASS: ',I,'out of ',NXP

+            CALL XABORT(HMSG)

+          ENDIF

+        ENDDO

+        LSPX=.TRUE.

+        GOTO 50

+      ELSEIF(TEXT.EQ.'SPLITY-ASS') THEN

+        DO I=1,NYP

+          CALL REDGET(INDIC,SYP(I),FLOT,TEXT,DFLOT)

+          IF(INDIC.NE.1) THEN

+            WRITE(HMSG,*) '@NAPGEO: integer number expected'

+     1  //' for SPLITY-ASS: ',I,'out of ',NYP

+            CALL XABORT(HMSG)

+          ENDIF

+        ENDDO

+        LSPY=.TRUE.

+        GOTO 50

+C read final ';'

+      ELSEIF(TEXT.EQ.'MAX-MIX-GEO') THEN

+        CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)

+        IF(INDIC.NE.1) CALL XABORT('@NAPGEO: integer data expected.')

+        NMIXD2=MAX(NMIXD,NITMA)

+        GOTO 50

+C read final ';'

+      ELSEIF(TEXT.EQ.';') THEN

+        if(impx.ge.5) then

+          write(6,*) 'Splitting within assemblies:'

+          write(6,*) 'SXP',(SXP(I),I=1,NXP)

+          write(6,*) 'SYP',(SYP(I),I=1,NYP)

+        endif

+        GOTO 60

+      ELSE

+        CALL XABORT('@NAPGEO: '//TEXT//' WRONG KEYWORD')

+      ENDIF

+C compute new dimension

+C     get number of x and y original mesh to split

+   60 NXA=0

+      DO I=1,NXD

+        NXA=NXA+AXD(I)

+      ENDDO

+      NYA=0

+      DO I=1,NYD

+        NYA=NYA+AYD(I)

+      ENDDO

+C     compute new dimension

+      NXF=NXD+NXA*(NXP-1)

+      NYF=NYD+NYA*(NYP-1)

+      NZF=NZD

+      NREGF=NXF*NYF*NZF

+C allocate new geometry dimensions

+      ALLOCATE(MXF(NXF+1),MYF(NYF+1))

+      ALLOCATE(MIXF(NXF,NYF,NZF))

+      ALLOCATE(SXF(NXF),SYF(NYF))

+C Compute new x/y mesh and new x/y split

+      if(impx.ge.100)write(6,*) 'debug: Compute new x/y mesh and split'

+      J=1

+      MXF(J)=MXD(1)

+      DO I=2,NXD+1

+        IF(AXD(I-1).EQ.0) THEN

+          J=J+1

+          MXF(J)=MXD(I)

+          SXF(J-1)=SXD(I-1)

+        ELSE

+          DO K=2,NXP+1

+            J=J+1

+            MXF(J)=MXD(I-1)+MXP(K)-MXP(1)

+            SXF(J-1)=SXP(K-1)

+          ENDDO

+        ENDIF

+      ENDDO

+      J=1

+      MYF(J)=MYD(1)

+      DO I=2,NYD+1

+        IF(AYD(I-1).EQ.0) THEN

+          J=J+1

+          MYF(J)=MYD(I)

+          SYF(J-1)=SYD(I-1)

+        ELSE

+          DO K=2,NYP+1

+            J=J+1

+            MYF(J)=MYD(I-1)+MYP(K)-MYP(1)

+            SYF(J-1)=SYP(K-1)

+          ENDDO

+        ENDIF

+      ENDDO

+      IF(MXF(NXF+1).NE.MXD(NXD+1)) CALL XABORT('@NAPGEO: OLD and NEW'

+     1 //' X MESH do not match.')

+      IF(MYF(NYF+1).NE.MYD(NYD+1)) CALL XABORT('@NAPGEO: OLD and NEW'

+     1 //' Y MESH do not match.')

+C Compute new mixture

+      if(impx.ge.100)write(6,*) 'debug: Compute new mixture'

+      NMIXF=NMIXD2+NMIXA*NMIXP

+      MIXF(:NXF,:NYF,:NZF)=-1

+      DO 100 L=1,NMIXA

+         MIXA(L+NMIXA)=NMIXD2+(L-1)*NMIXP+1 

+  100 CONTINUE

+      DO 240 K=1,NZD

+      JFBEG=0

+      DO 230 J=1,NYD

+      IFBEG=0

+      DO 220 I=1,NXD

+        LMIX=0

+        DO 110 L=1,NMIXA

+          IF(MIXA(L).EQ.MIXD(I,J,K)) THEN

+            LMIX=L

+          ENDIF

+  110   CONTINUE

+        IF((AXD(I).EQ.1).AND.(AYD(J).EQ.1).AND.(LMIX.NE.0)) THEN

+          DO 130 JP=1,NYP

+          JF=JFBEG+JP

+          DO 120 IP=1,NXP

+            MIXF(IFBEG+IP,JF,K)=MIXA(LMIX+NMIXA)-1+MIXP(IP,JP)

+  120     CONTINUE

+  130     CONTINUE

+          IFBEG=IFBEG+NXP

+        ELSE

+          NXPTMP=1

+          IF(AXD(I).EQ.1) NXPTMP=NXP

+          NYPTMP=1

+          IF(AYD(J).EQ.1) NYPTMP=NYP

+          DO 150 JP=1,NYPTMP

+          JF=JFBEG+JP

+          DO 140 IP=1,NXPTMP

+            MIXF(IFBEG+IP,JF,K)=MIXD(I,J,K)

+  140     CONTINUE

+  150     CONTINUE

+          IFBEG=IFBEG+NXPTMP

+        ENDIF

+  220 CONTINUE

+      NYPTMP=1

+      IF(AYD(J).EQ.1) NYPTMP=NYP

+      JFBEG=JFBEG+NYPTMP

+  230 CONTINUE

+  240 CONTINUE

+      

+C Compute A-ZONE

+      if(impx.ge.100)write(6,*) 'debug: Compute A-ZONE'

+      IASS=0

+      ALLOCATE(NBAX(NYD))

+      ALLOCATE(IBAX(NYD))

+      DO 340 J=1,NYD

+      NBAX(J)=0

+      IBAX(J)=0

+      DO 330 I=1,NXD

+      DO 320 K=1,NZD

+        DO 310 L=1,NMIXA

+          IF(MIXA(L).EQ.MIXD(I,J,K)) THEN

+            IASS=IASS+1

+            NBAX(J)=NBAX(J)+1

+            IF(IBAX(J).EQ.0) IBAX(J)=I

+            GOTO 330

+          ENDIF

+  310   CONTINUE

+  320 CONTINUE

+  330 CONTINUE

+  340 CONTINUE

+      NASS=IASS

+*

+      ALLOCATE(AZONE(NASS*NXP*NYP))

+      IZ=0

+      IASS=0

+      DO 370 J=1,NYD

+      DO 360 I=1,NBAX(J)

+      IASS=IASS+1

+      DO 365 JP=1,NYP

+      DO 355 IP=1,NXP

+      IZT=IZ+(JP-1)*NXP*NBAX(J)+(I-1)*NXP+IP

+      AZONE(IZT)=IASS

+  355 CONTINUE

+  365 CONTINUE

+  360 CONTINUE

+      IZ=IZ+NBAX(J)*NXP*NYP

+  370 CONTINUE

+*      

+      if(impx.ge.5)then

+        write(6,*) 'New mixtures:'

+        do K=1,NZF

+          write(6,*) 'plane #',K

+          do J=1,NYF

+            write(6,*) (MIXF(I,J,K),I=1,NXF)

+          enddo

+        enddo

+

+        write(6,*) 'Assembly zones:'

+        IZ=0

+        do J=1,NYD

+          do K=1,NYP

+            write(6,*) (AZONE(I),I=IZ+(K-1)*NBAX(J)*NXP+1,

+     1      IZ+K*NBAX(J)*NXP)

+          enddo

+          IZ=IZ+NBAX(J)*NXP*NYP

+        enddo

+      endif

+C Verify new mixture

+      DO K=1,NZF

+      DO J=1,NYF

+      DO I=1,NXF

+      IF(MIXF(I,J,K).EQ.-1) CALL XABORT('@NAPGEO: new '

+     1 //'geometry mixture not assigned')

+      ENDDO

+      ENDDO

+      ENDDO

+

+      DEALLOCATE(MXD,MYD)

+      DEALLOCATE(SXD,SYD)

+      DEALLOCATE(MIXD)

+

+      DEALLOCATE(MXP,MYP)

+      DEALLOCATE(SXP,SYP)

+      DEALLOCATE(MIXP)

+C Compute relative position of assembly in original geometry

+      if(impx.ge.100)write(6,*) 'debug: Compute relative position'

+      JM=0

+      JN=0

+      LPOS=.TRUE.

+      DO J=1,NYD

+        IF(NBAX(J).NE.0) THEN

+          JN=JN+1

+          IF(LPOS) THEN

+            JM=J

+            LPOS=.FALSE.

+          ENDIF

+        ENDIF

+      ENDDO

+

+      IM=10000000

+      DO J=JM,JM+JN-1

+        IM=MIN(IBAX(J),IM)

+      ENDDO

+      DO J=1,NYD

+        IBAX(J)=IBAX(J)-IM+1

+      ENDDO

+

+C Save heterogeneous core geometry

+      if(impx.ge.100)write(6,*) 'debug: Save heter. core geometry'

+      ISTATE(:NSTATE)=0

+      ISTATE(1)=ITYPGP

+      ISTATE(3)=NXF

+      ISTATE(4)=NYF

+      ISTATE(5)=NZF

+      ISTATE(6)=NREGF

+      ISTATE(7)=NMIXF

+*      ISTATE(39)=NMIXA

+*      ISTATE(40)=NMIXP

+      IF(LSPX .OR. LSPY .OR. LSPZ) ISTATE(11)=1

+      CALL LCMPUT(IPGNW,'STATE-VECTOR',NSTATE,1,ISTATE)

+      CALL LCMPUT(IPGNW,'MESHX',NXF+1,2,MXF)

+      CALL LCMPUT(IPGNW,'MESHY',NYF+1,2,MYF)

+      CALL LCMPUT(IPGNW,'MESHZ',NZF+1,2,MZD)

+      IF(LSPX) CALL LCMPUT(IPGNW,'SPLITX',NXF,1,SXF)

+      IF(LSPY) CALL LCMPUT(IPGNW,'SPLITY',NYF,1,SYF)

+      IF(LSPZ) CALL LCMPUT(IPGNW,'SPLITZ',NZF,1,SZD)

+      CALL LCMPUT(IPGNW,'MIX',NREGF,1,MIXF)

+      CALL LCMPUT(IPGNW,'NCODE',6,1,NCODE)

+      CALL LCMPUT(IPGNW,'ICODE',6,1,ICODE)

+      CALL LCMPUT(IPGNW,'ZCODE',6,2,ZCODE)

+      CALL LCMPUT(IPGNW,'MIX-ASBLY',2*NMIXA,1,MIXA)

+      CALL LCMPUT(IPGNW,'SIGNATURE',3,3,KCHAR)

+      CALL LCMPUT(IPGNW,'A-ZONE',NASS*NXP*NYP,1,AZONE)

+      CALL LCMPUT(IPGNW,'A-NX',JN,1,NBAX(JM))

+      CALL LCMPUT(IPGNW,'A-IBX',JN,1,IBAX(JM))

+      CALL LCMPUT(IPGNW,'A-NMIXP',1,1,NMIXP)

+! 

+      if(impx.ge.100)write(6,*) 'debug: beging deallacate'

+

+      DEALLOCATE(IBAX)

+      DEALLOCATE(NBAX)

+      DEALLOCATE(AZONE)

+      DEALLOCATE(MXF,MYF)

+      DEALLOCATE(MZD)

+      DEALLOCATE(SXF,SYF)

+      DEALLOCATE(SZD)

+      DEALLOCATE(MIXF)

+      DEALLOCATE(MIXA)

+

+

+      RETURN

+      END