Initial commit from Polytechnique Montreal

1 files changed, 145 insertions, 0 deletions

diff --git a/Trivac/src/FLDORT.f b/Trivac/src/FLDORT.f
new file mode 100755
index 0000000..41d7ff0
--- /dev/null
+++ b/Trivac/src/FLDORT.f
@@ -0,0 +1,145 @@
+*DECK FLDORT
+      SUBROUTINE FLDORT(IPSYS,IPFLUX,NUN,NGRP,LMOD)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Test the biorthogonality of the direct-CADjoint eigenvectors.
+*
+*Copyright:
+* Copyright (C) 2020 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
+* IPSYS   L_SYSTEM pointer to system matrices.
+* IPFLUX  L_FLUX pointer to the solution.
+* NUN     number of unknowns in each energy group.
+* NGRP    number of energy groups.
+* LMOD    number of orthogonal harmonics to compute.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPSYS,IPFLUX
+      INTEGER NUN,NGRP,LMOD
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(NSTATE=40)
+      INTEGER ISTATE(NSTATE)
+      CHARACTER TEXT12*12,HSMG*131
+      TYPE(C_PTR) JPFLUX,KPFLUX,MPFLUX
+      REAL, DIMENSION(:), POINTER :: AGARM
+      TYPE(C_PTR) AGARM_PTR
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, DIMENSION(:), ALLOCATABLE :: GAR
+      COMPLEX, DIMENSION(:,:,:), ALLOCATABLE :: CEV,CAD
+      COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: DWORK,ORTHO
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(DWORK(NUN,NGRP),CEV(NUN,NGRP,LMOD),CAD(NUN,NGRP,LMOD),
+     1 ORTHO(LMOD,LMOD),GAR(NUN))
+*----
+*  FLUX RECOVERY
+*----
+      CALL LCMLEN(IPFLUX,'MODE',ILONG,ITYLCM)
+      IF((ILONG.EQ.0).AND.(LMOD.EQ.1)) THEN
+        MPFLUX=LCMGID(IPFLUX,'AFLUX')
+        DO IGR=1,NGRP
+          CALL LCMGDL(MPFLUX,IGR,GAR)
+          CAD(:NUN,IGR,1)=GAR(:NUN)
+        ENDDO
+        MPFLUX=LCMGID(IPFLUX,'FLUX')
+        DO IGR=1,NGRP
+          CALL LCMGDL(MPFLUX,IGR,GAR)
+          CEV(:NUN,IGR,1)=GAR(:NUN)
+        ENDDO
+      ELSE IF(ILONG.GT.0) THEN
+        DO IMOD=1,LMOD
+          JPFLUX=LCMGID(IPFLUX,'MODE')
+          CALL LCMLEL(JPFLUX,IMOD,ILONG,ITYLCM)
+          IF(ILONG.EQ.0) THEN
+            WRITE(6,'(20HFLDORT: MISSING MODE,I4,1H.)') IMOD
+            CALL XABORT(HSMG)
+          ENDIF
+          KPFLUX=LCMGIL(JPFLUX,IMOD)
+          MPFLUX=LCMGID(KPFLUX,'AFLUX')
+          DO IGR=1,NGRP
+            CALL LCMLEL(MPFLUX,IGR,ILONG,ITYLCM)
+            IF(ITYLCM.EQ.2) THEN
+              CALL LCMGDL(MPFLUX,IGR,GAR)
+              CAD(:NUN,IGR,IMOD)=GAR(:NUN)
+            ELSE IF(ITYLCM.EQ.6) THEN
+              CALL LCMGDL(MPFLUX,IGR,CAD(1,IGR,IMOD))
+            ENDIF
+          ENDDO
+          MPFLUX=LCMGID(KPFLUX,'FLUX')
+          DO IGR=1,NGRP
+            CALL LCMLEL(MPFLUX,IGR,ILONG,ITYLCM)
+            IF(ITYLCM.EQ.2) THEN
+              CALL LCMGDL(MPFLUX,IGR,GAR)
+              CEV(:NUN,IGR,IMOD)=GAR(:NUN)
+            ELSE IF(ITYLCM.EQ.6) THEN
+              CALL LCMGDL(MPFLUX,IGR,CEV(1,IGR,IMOD))
+            ENDIF
+          ENDDO
+        ENDDO
+      ELSE
+        CALL XABORT('FLDORT: MODE INFORMATION MISSING.')
+      ENDIF
+*----
+*  MULTIPLY FLUX WITH B MATRIX
+*----
+      CALL LCMGET(IPSYS,'STATE-VECTOR',ISTATE)
+      LL4=ISTATE(2)
+      DO JMOD=1,LMOD
+        DWORK(:NUN,:NGRP)=0.0D0
+        DO IGR=1,NGRP
+          DO JGR=1,NGRP
+            WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+            CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+            IF(ILONG.EQ.0) CYCLE
+            CALL LCMGPD(IPSYS,TEXT12,AGARM_PTR)
+            CALL C_F_POINTER(AGARM_PTR,AGARM,(/ ILONG /))
+            DO I=1,ILONG
+              DWORK(I,IGR)=DWORK(I,IGR)+CMPLX(AGARM(I)*CEV(I,JGR,JMOD))
+            ENDDO
+          ENDDO
+        ENDDO
+*----
+*  COMPUTE ORTHONORMAL MATRIX
+*----
+        DO IMOD=1,LMOD
+          ORTHO(IMOD,JMOD)=0.0D0
+          DO I=1,LL4
+            DO IGR=1,NGRP
+              ORTHO(IMOD,JMOD)=ORTHO(IMOD,JMOD)+CAD(I,IGR,IMOD)*
+     1        DWORK(I,IGR)
+            ENDDO
+          ENDDO
+        ENDDO
+      ENDDO
+*----
+*  PRINT ORTHONORMAL MATRIX
+*----
+      WRITE(6,'(/28H FLDORT: ORTHONORMAL MATRIX:)')
+      DO IMOD=1,LMOD
+        WRITE(6,'(3X,1P,15E12.4)') REAL(ORTHO(IMOD,:LMOD))
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(GAR,ORTHO,CAD,CEV,DWORK)
+      RETURN
+      END