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 --- Trivac/src/GPTGRA.f | 298 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 298 insertions(+) create mode 100755 Trivac/src/GPTGRA.f (limited to 'Trivac/src/GPTGRA.f') 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 -- cgit v1.2.3