summaryrefslogtreecommitdiff
path: root/Trivac/src/GPTGRA.f
diff options
context:
space:
mode:
Diffstat (limited to 'Trivac/src/GPTGRA.f')
-rwxr-xr-xTrivac/src/GPTGRA.f298
1 files changed, 298 insertions, 0 deletions
diff --git a/Trivac/src/GPTGRA.f b/Trivac/src/GPTGRA.f
new file mode 100755
index 0000000..f1867a3
--- /dev/null
+++ b/Trivac/src/GPTGRA.f
@@ -0,0 +1,298 @@
+*DECK GPTGRA
+ SUBROUTINE GPTGRA(IPTRK,IPSYS,IPFLUP,LADJ,LGAR1,LL4,ITY,NUN,NGRP,
+ 1 ICL1,ICL2,IMPX,NNADI,MAXINR,EPSINR,EVAL,EVECT,ADECT,EASS,SOUR,
+ 2 GAR1,ITER,TKT,TKB,ZNORM,GRAD)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute multigroup delta flux in a fixed source eigenvalue iteration.
+*
+*Copyright:
+* Copyright (C) 2019 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
+* IPTRK L_TRACK pointer to the tracking information.
+* IPSYS L_SYSTEM pointer to system matrices.
+* IPFLUP L_FLUX pointer to the gpt solution
+* LADJ flag set to .TRUE. for adjoint solution acceleration.
+* LGAR1 flag set to .TRUE. for recomputing GAR1.
+* LL4 order of the system matrices.
+* ITY type of solution (2: classical Trivac; 3: Thomas-Raviart).
+* NUN number of unknowns in each energy group.
+* NGRP number of energy groups.
+* ICL1 number of free up-scattering iterations in one cycle of the
+* inverse power method.
+* ICL2 number of accelerated up-scattering iterations in one cycle.
+* IMPX print parameter (set to 0 for no printing).
+* NNADI number of inner ADI iterations per outer iteration.
+* MAXINR maximum number of thermal iterations.
+* EPSINR thermal iteration epsilon.
+* EVAL eigenvalue.
+* EVECT unknown vector for the non perturbed direct flux
+* ADECT unknown vector for the non perturbed adjoint flux
+* EASS solution of the fixed source eigenvalue problem
+* SOUR fixed source
+* GAR1 delta flux for this iteration before Hotelling deflation.
+*
+*Parameters: input/output
+* ITER actual number of thermal iterations.
+* TKT CPU time spent to compute the solution of linear systems.
+* TKB CPU time spent to compute the bilinear products.
+* ZNORM Hotelling deflation accuracy.
+* GRAD delta flux for this iteration.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK,IPSYS,IPFLUP
+ LOGICAL LADJ,LGAR1
+ INTEGER LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,NNADI,MAXINR,ITER
+ REAL EPSINR,EVECT(NUN,NGRP),ADECT(NUN,NGRP),EASS(NUN,NGRP),
+ 1 SOUR(NUN,NGRP),GAR1(NUN,NGRP),TKT,TKB,GRAD(NUN,NGRP)
+ DOUBLE PRECISION EVAL,ZNORM
+*----
+* LOCAL VARIABLES
+*----
+ CHARACTER*12 TEXT12
+ DOUBLE PRECISION DDELN1,DDELD1
+ REAL, DIMENSION(:), ALLOCATABLE :: WORK1,WORK3
+ REAL, DIMENSION(:), POINTER :: AGAR
+ TYPE(C_PTR) AGAR_PTR
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(WORK1(NUN))
+*
+ IF(LADJ) THEN
+ CALL KDRCPU(TK1)
+* ADJOINT SOLUTION
+ IF(LGAR1) THEN
+ DO 55 IGR=1,NGRP
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EASS(1,IGR),
+ 1 GAR1(1,IGR))
+ DO 50 JGR=1,NGRP
+ IF(JGR.EQ.IGR) GO TO 30
+ WRITE(TEXT12,'(1HA,2I3.3)') JGR,IGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 30
+ IF(ITY.EQ.13) THEN
+ ALLOCATE(WORK3(LL4))
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EASS(1,JGR),
+ 1 WORK3(1))
+ DO 10 I=1,LL4
+ GAR1(I,IGR)=GAR1(I,IGR)-WORK3(I)
+ 10 CONTINUE
+ DEALLOCATE(WORK3)
+ ELSE
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 20 I=1,ILONG
+ GAR1(I,IGR)=GAR1(I,IGR)-AGAR(I)*EASS(I,JGR)
+ 20 CONTINUE
+ ENDIF
+ 30 WRITE(TEXT12,'(1HB,2I3.3)') JGR,IGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 50
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 40 I=1,ILONG
+ GAR1(I,IGR)=GAR1(I,IGR)-REAL(EVAL)*AGAR(I)*EASS(I,JGR)
+ 40 CONTINUE
+ 50 CONTINUE
+ 55 CONTINUE
+ ENDIF
+*----
+* DIRECTION EVALUATION.
+*----
+ DO 100 IGR=NGRP,1,-1
+ DO 60 I=1,LL4
+ GRAD(I,IGR)=-SOUR(I,IGR)-GAR1(I,IGR)
+ 60 CONTINUE
+ DO 90 JGR=NGRP,IGR+1,-1
+ WRITE(TEXT12,'(1HA,2I3.3)') JGR,IGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 90
+ IF(ITY.EQ.13) THEN
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,JGR),WORK1(1))
+ DO 70 I=1,LL4
+ GRAD(I,IGR)=GRAD(I,IGR)+WORK1(I)
+ 70 CONTINUE
+ ELSE
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 80 I=1,ILONG
+ GRAD(I,IGR)=GRAD(I,IGR)+AGAR(I)*GRAD(I,JGR)
+ 80 CONTINUE
+ ENDIF
+ 90 CONTINUE
+*
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+ CALL FLDADI(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,IGR),NNADI)
+ 100 CONTINUE
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+*----
+* PERFORM THERMAL (UP-SCATTERING) ITERATIONS
+*----
+ ITER=1
+ IF(MAXINR.GT.1) THEN
+ CALL FLDTHR(IPTRK,IPSYS,IPFLUP,.TRUE.,LL4,ITY,NUN,NGRP,
+ 1 ICL1,ICL2,IMPX,NNADI,0,MAXINR,EPSINR,ITER,TKT,TKB,GRAD)
+ ENDIF
+*----
+* HOTELLING DEFLATION.
+*----
+ CALL KDRCPU(TK1)
+ DDELN1=0.0D0
+ DDELD1=0.0D0
+ DO 135 IGR=1,NGRP
+ WORK1(:LL4)=0.0
+ DO 120 JGR=1,NGRP
+ WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 120
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 110 I=1,ILONG
+ WORK1(I)=WORK1(I)+AGAR(I)*EVECT(I,JGR)
+ 110 CONTINUE
+ 120 CONTINUE
+ DO 130 I=1,LL4
+ DDELN1=DDELN1+WORK1(I)*EASS(I,IGR)
+ DDELD1=DDELD1+WORK1(I)*ADECT(I,IGR)
+ 130 CONTINUE
+ 135 CONTINUE
+ ZNORM=DDELN1/DDELD1
+ DO 145 IGR=1,NGRP
+ DO 140 I=1,LL4
+ GRAD(I,IGR)=GRAD(I,IGR)-REAL(ZNORM)*ADECT(I,IGR)
+ 140 CONTINUE
+ 145 CONTINUE
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+ ELSE
+ CALL KDRCPU(TK1)
+* DIRECT SOLUTION
+ IF(LGAR1) THEN
+ DO 195 IGR=1,NGRP
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EASS(1,IGR),
+ 1 GAR1(1,IGR))
+ DO 190 JGR=1,NGRP
+ IF(JGR.EQ.IGR) GO TO 170
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 170
+ IF(ITY.EQ.13) THEN
+ ALLOCATE(WORK3(LL4))
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EASS(1,JGR),
+ 1 WORK3(1))
+ DO 150 I=1,LL4
+ GAR1(I,IGR)=GAR1(I,IGR)-WORK3(I)
+ 150 CONTINUE
+ DEALLOCATE(WORK3)
+ ELSE
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 160 I=1,ILONG
+ GAR1(I,IGR)=GAR1(I,IGR)-AGAR(I)*EASS(I,JGR)
+ 160 CONTINUE
+ ENDIF
+ 170 WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 190
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 180 I=1,ILONG
+ GAR1(I,IGR)=GAR1(I,IGR)-REAL(EVAL)*AGAR(I)*EASS(I,JGR)
+ 180 CONTINUE
+ 190 CONTINUE
+ 195 CONTINUE
+ ENDIF
+*----
+* DIRECTION EVALUATION.
+*----
+ DO 240 IGR=1,NGRP
+ DO 200 I=1,LL4
+ GRAD(I,IGR)=-SOUR(I,IGR)-GAR1(I,IGR)
+ 200 CONTINUE
+ DO 230 JGR=1,IGR-1
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 230
+ IF(ITY.EQ.13) THEN
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,JGR),WORK1(1))
+ DO 210 I=1,LL4
+ GRAD(I,IGR)=GRAD(I,IGR)+WORK1(I)
+ 210 CONTINUE
+ ELSE
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 220 I=1,ILONG
+ GRAD(I,IGR)=GRAD(I,IGR)+AGAR(I)*GRAD(I,JGR)
+ 220 CONTINUE
+ ENDIF
+ 230 CONTINUE
+*
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+ CALL FLDADI(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,IGR),NNADI)
+ 240 CONTINUE
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+*----
+* PERFORM THERMAL (UP-SCATTERING) ITERATIONS
+*----
+ ITER=1
+ IF(MAXINR.GT.1) THEN
+ CALL FLDTHR(IPTRK,IPSYS,IPFLUP,.FALSE.,LL4,ITY,NUN,NGRP,
+ 1 ICL1,ICL2,IMPX,NNADI,0,MAXINR,EPSINR,ITER,TKT,TKB,GRAD)
+ ENDIF
+*----
+* HOTELLING DEFLATION.
+*----
+ CALL KDRCPU(TK1)
+ DDELN1=0.0D0
+ DDELD1=0.0D0
+ DO 275 IGR=1,NGRP
+ WORK1(:LL4)=0.0
+ DO 260 JGR=1,NGRP
+ WRITE(TEXT12,'(1HB,2I3.3)') JGR,IGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 260
+ CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+ CALL C_F_POINTER(AGAR_PTR,AGAR,(/ NUN /))
+ DO 250 I=1,ILONG
+ WORK1(I)=WORK1(I)+AGAR(I)*ADECT(I,JGR)
+ 250 CONTINUE
+ 260 CONTINUE
+ DO 270 I=1,LL4
+ DDELN1=DDELN1+WORK1(I)*EASS(I,IGR)
+ DDELD1=DDELD1+WORK1(I)*EVECT(I,IGR)
+ 270 CONTINUE
+ 275 CONTINUE
+ ZNORM=DDELN1/DDELD1
+ DO 285 IGR=1,NGRP
+ DO 280 I=1,LL4
+ GRAD(I,IGR)=GRAD(I,IGR)-REAL(ZNORM)*EVECT(I,IGR)
+ 280 CONTINUE
+ 285 CONTINUE
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+ ENDIF
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(WORK1)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 15:15:38 +0000