Initial commit from Polytechnique Montreal

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