summaryrefslogtreecommitdiff
path: root/Dragon/src/PSOISO.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/PSOISO.f')
-rw-r--r--Dragon/src/PSOISO.f295
1 files changed, 295 insertions, 0 deletions
diff --git a/Dragon/src/PSOISO.f b/Dragon/src/PSOISO.f
new file mode 100644
index 0000000..9782bcf
--- /dev/null
+++ b/Dragon/src/PSOISO.f
@@ -0,0 +1,295 @@
+*DECK PSOISO
+ SUBROUTINE PSOISO(IPTRK,IPGEOM,NREG,LX,LY,LZ,NG,NUNS,NDIM,
+ 1 NSOUR,ISOUR,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,XXX,YYY,ZZZ,MESHL,
+ 2 SUNKNO,NORM)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute moments of fixed isotropic sources.
+*
+*Copyright:
+* Copyright (C) 2022 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): C. Bienvenue
+*
+*Parameters: input
+* IPTRK pointer to the tracking LCM object.
+* IPGEOM pointer to the geometry LCM object.
+* NREG number of regions.
+* LX number of meshes along X axis.
+* LY number of meshes along Y axis.
+* LZ number of meshes along Z axis.
+* NG number of energy groups.
+* NUNS number of unknowns in vector SUNKNO.
+* NDIM geometry dimension.
+* NSOUR number of sources defined.
+* ISOUR intensity of the sources.
+* XMIN lower boundaries of the sources along X axis.
+* XMAX upper boundaries of the sources along X axis.
+* YMIN lower boundaries of the sources along Y axis.
+* YMAX upper boundaries of the sources along Y axis.
+* ZMIN lower boundaries of the sources along Z axis.
+* ZMAX upper boundaries of the sources along Z axis.
+* XXX regions boundaries along X axis.
+* YYY regions boundaries along Y axis.
+* ZZZ regions boundaries along Z axis.
+* MESHL number of regions along X-, Y- and Z-axis
+*
+*Parameters: output
+* SUNKNO source vector.
+* NORM normalization factor.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK,IPGEOM
+ INTEGER NREG,LX,LY,LZ,NG,NUNS,NDIM,NSOUR,MESHL(3)
+ REAL XMIN(NSOUR),XMAX(NSOUR),YMIN(NSOUR),YMAX(NSOUR),ZMIN(NSOUR),
+ 1 ZMAX(NSOUR),ISOUR(NG),SUNKNO(NUNS,NG),XXX(MESHL(1)),
+ 2 YYY(MESHL(2)),ZZZ(MESHL(3)),NORM
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER(NSTATE=40)
+ INTEGER ISTATE(NSTATE),SPLIT_LEN(3),XP(NREG),YP(NREG),ZP(NREG),
+ 1 EELEM
+ REAL X(LX),Y(LY),Z(LZ)
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER,ALLOCATABLE,DIMENSION(:) :: SPLITX,SPLITY,SPLITZ
+*----
+* RECOVER TRACKING INFORMATION
+*----
+ CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+ ITYPE=ISTATE(6)
+ NSCT=ISTATE(7)
+ IELEM=ISTATE(8)
+ ISCAT=ISTATE(16)
+ EELEM=ISTATE(35)
+*----
+* RECOVER GEOMETRY INFORMATION
+*----
+ IF(NDIM.EQ.1) THEN
+ CALL LCMLEN(IPGEOM,'SPLITX',SPLIT_LEN(1),ITYLCM)
+ ALLOCATE(SPLITX(SPLIT_LEN(1)))
+ CALL LCMGET(IPGEOM,'SPLITX',SPLITX)
+ ELSE IF(NDIM.EQ.2) THEN
+ CALL LCMLEN(IPGEOM,'SPLITX',SPLIT_LEN(1),ITYLCM)
+ CALL LCMLEN(IPGEOM,'SPLITY',SPLIT_LEN(2),ITYLCM)
+ ALLOCATE(SPLITX(SPLIT_LEN(1)),SPLITY(SPLIT_LEN(2)))
+ CALL LCMGET(IPGEOM,'SPLITX',SPLITX)
+ CALL LCMGET(IPGEOM,'SPLITY',SPLITY)
+ ELSE IF(NDIM.EQ.3) THEN
+ CALL LCMLEN(IPGEOM,'SPLITX',SPLIT_LEN(1),ITYLCM)
+ CALL LCMLEN(IPGEOM,'SPLITY',SPLIT_LEN(2),ITYLCM)
+ CALL LCMLEN(IPGEOM,'SPLITZ',SPLIT_LEN(3),ITYLCM)
+ ALLOCATE(SPLITX(SPLIT_LEN(1)),SPLITY(SPLIT_LEN(2)),
+ 1 SPLITZ(SPLIT_LEN(3)))
+ CALL LCMGET(IPGEOM,'SPLITX',SPLITX)
+ CALL LCMGET(IPGEOM,'SPLITY',SPLITY)
+ CALL LCMGET(IPGEOM,'SPLITZ',SPLITZ)
+ ENDIF
+*----
+* 1D CARTESIAN CASE
+*----
+
+ IF(NDIM.EQ.1) THEN
+
+ ! CALCULATE X-COORDINATES OF EACH VOXELS
+ K=1
+ DO I=1,SPLIT_LEN(1)
+ DO J=1,SPLITX(I)
+ XP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+
+ DO 10 IX=1,LX
+ STEPX=(XXX(XP(IX)+1)-XXX(XP(IX)))/SPLITX(XP(IX))
+ IF(XP(IX).EQ.1) THEN
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-1)
+ ELSE
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-SUM(SPLITX(1:XP(IX)-1))-1)
+ ENDIF
+ 10 CONTINUE
+
+ ! CALCULATE THE SOURCE DENSITY
+ NORM=0.0
+ DO 40 IX=1,LX
+ IR=IX
+ DO 30 N=1,NSOUR
+ IF(XMIN(N).LE.X(IX).AND.XMAX(N).GE.X(IX)) THEN
+ IND=(IR-1)*NSCT*IELEM*EELEM+1
+ DO 20 IG=1,NG
+ SUNKNO(IND,IG)=SUNKNO(IND,IG)+ISOUR(IG)
+ IF(N.EQ.NSOUR) THEN
+ NORM=NORM+(XXX(XP(IX)+1)-XXX(XP(IX)))/SPLITX(XP(IX))*
+ 1 SUNKNO(IND,IG)
+ ENDIF
+ 20 CONTINUE
+ ENDIF
+ 30 CONTINUE
+ 40 CONTINUE
+
+*----
+* 2D CARTESIAN CASE
+*----
+
+ ELSE IF(NDIM.EQ.2) THEN
+
+ ! CALCULATE XY-COORDINATES OF EACH VOXELS
+ K=1
+ DO I=1,SPLIT_LEN(1)
+ DO J=1,SPLITX(I)
+ XP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+ K=1
+ DO I=1,SPLIT_LEN(2)
+ DO J=1,SPLITY(I)
+ YP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+
+ DO 100 IX=1,LX
+ STEPX=(XXX(XP(IX)+1)-XXX(XP(IX)))/SPLITX(XP(IX))
+ IF(XP(IX).EQ.1) THEN
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-1)
+ ELSE
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-SUM(SPLITX(1:XP(IX)-1))-1)
+ ENDIF
+ 100 CONTINUE
+
+ DO 110 IY=1,LY
+ STEPY=(YYY(YP(IY)+1)-YYY(YP(IY)))/SPLITY(YP(IY))
+ IF(YP(IY).EQ.1) THEN
+ Y(IY)=YYY(YP(IY))+0.5*STEPY+STEPY*(IY-1)
+ ELSE
+ Y(IY)=YYY(YP(IY))+0.5*STEPY+STEPY*(IY-SUM(SPLITY(1:YP(IY)-1))-1)
+ ENDIF
+ 110 CONTINUE
+
+ ! CALCULATE THE SOURCE DENSITY
+ NORM=0.0
+ DO 150 IY=1,LY
+ DO 140 IX=1,LX
+ IR=IX+(IY-1)*LX
+ DO 130 N=1,NSOUR
+ IF(XMIN(N).LE.X(IX).AND.XMAX(N).GE.X(IX).AND.
+ 1 YMIN(N).LE.Y(IY).AND.YMAX(N).GE.Y(IY)) THEN
+ IND=(IR-1)*NSCT*IELEM*IELEM*EELEM+1
+ DO 120 IG=1,NG
+ SUNKNO(IND,IG)=SUNKNO(IND,IG)+ISOUR(IG)
+ IF(N.EQ.NSOUR) THEN
+ NORM=NORM+( (XXX(XP(IX)+1)-XXX(XP(IX)))/SPLITX(XP(IX))*
+ 1 (YYY(YP(IY)+1)-YYY(YP(IY)))/SPLITY(YP(IY)))*SUNKNO(IND,IG)
+ ENDIF
+ 120 CONTINUE
+ ENDIF
+ 130 CONTINUE
+ 140 CONTINUE
+ 150 CONTINUE
+
+*----
+* 3D CARTESIAN CASE
+*----
+
+ ELSE IF(NDIM.EQ.3) THEN
+
+ ! CALCULATE XYZ-COORDINATES OF EACH VOXELS
+ K=1
+ DO I=1,SPLIT_LEN(1)
+ DO J=1,SPLITX(I)
+ XP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+ K=1
+ DO I=1,SPLIT_LEN(2)
+ DO J=1,SPLITY(I)
+ YP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+ K=1
+ DO I=1,SPLIT_LEN(3)
+ DO J=1,SPLITZ(I)
+ ZP(K)=I
+ K=K+1
+ ENDDO
+ ENDDO
+
+ DO 200 IX=1,LX
+ STEPX=(XXX(XP(IX)+1)-XXX(XP(IX)))/SPLITX(XP(IX))
+ IF(XP(IX).EQ.1) THEN
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-1)
+ ELSE
+ X(IX)=XXX(XP(IX))+0.5*STEPX+STEPX*(IX-SUM(SPLITX(1:XP(IX)-1))-1)
+ ENDIF
+ 200 CONTINUE
+
+ DO 210 IY=1,LY
+ STEPY=(YYY(YP(IY)+1)-YYY(YP(IY)))/SPLITY(YP(IY))
+ IF(YP(IY).EQ.1) THEN
+ Y(IY)=YYY(YP(IY))+0.5*STEPY+STEPY*(IY-1)
+ ELSE
+ Y(IY)=YYY(YP(IY))+0.5*STEPY+STEPY*(IY-SUM(SPLITY(1:YP(IY)-1))-1)
+ ENDIF
+ 210 CONTINUE
+
+ DO 220 IZ=1,LZ
+ STEPZ=(ZZZ(ZP(IZ)+1)-ZZZ(ZP(IZ)))/SPLITZ(ZP(IZ))
+ IF(ZP(IZ).EQ.1) THEN
+ Z(IZ)=ZZZ(ZP(IZ))+0.5*STEPZ+STEPZ*(IZ-1)
+ ELSE
+ Z(IZ)=ZZZ(ZP(IZ))+0.5*STEPZ+STEPZ*(IZ-SUM(SPLITZ(1:ZP(IZ)-1))-1)
+ ENDIF
+ 220 CONTINUE
+
+ ! CALCULATE THE SOURCE DENSITY
+ NORM=0.0
+ DO 270 IZ=1,LZ
+ DO 260 IY=1,LY
+ DO 250 IX=1,LX
+ IR=IX+(IY-1)*LX+(IZ-1)*LX*LY
+ DO 240 N=1,NSOUR
+ IF(XMIN(N).LE.X(IX).AND.XMAX(N).GE.X(IX).AND.
+ 1 YMIN(N).LE.Y(IY).AND.YMAX(N).GE.Y(IY).AND.
+ 2 ZMIN(N).LE.Z(IZ).AND.ZMAX(N).GE.Z(IZ)) THEN
+ IND=(IR-1)*NSCT*IELEM*IELEM*IELEM*EELEM+1
+ DO 230 IG=1,NG
+ SUNKNO(IND,IG)=SUNKNO(IND,IG)+ISOUR(IG)
+ IF(N.EQ.NSOUR) THEN
+ NORM=NORM+((XXX(XP(IX)+1)-XXX(XP(IX)))/SPLITX(XP(IX))*
+ 1 (YYY(YP(IY)+1)-YYY(YP(IY)))/SPLITY(YP(IY))*
+ 2 (ZZZ(ZP(IZ)+1)-ZZZ(ZP(IZ)))/SPLITZ(ZP(IZ)))*SUNKNO(IND,IG)
+ ENDIF
+ 230 CONTINUE
+ ENDIF
+ 240 CONTINUE
+ 250 CONTINUE
+ 260 CONTINUE
+ 270 CONTINUE
+
+ ELSE
+ CALL XABORT('SOUR: INVALID GEOMETRY, ONLY 1D, 2D AND 3D CARTESIAN'
+ 1 //' GEOMETRY ARE ACTUALLY IMPLEMENTED.')
+ ENDIF
+
+ IF(ALLOCATED(SPLITX)) DEALLOCATE(SPLITX)
+ IF(ALLOCATED(SPLITY)) DEALLOCATE(SPLITY)
+ IF(ALLOCATED(SPLITZ)) DEALLOCATE(SPLITZ)
+
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:33:40 +0000