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 --- Donjon/src/NAPGEO.f | 487 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 487 insertions(+) create mode 100644 Donjon/src/NAPGEO.f (limited to 'Donjon/src/NAPGEO.f') 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 -- cgit v1.2.3