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/FLDTHR.f | 300 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 300 insertions(+) create mode 100755 Trivac/src/FLDTHR.f (limited to 'Trivac/src/FLDTHR.f') diff --git a/Trivac/src/FLDTHR.f b/Trivac/src/FLDTHR.f new file mode 100755 index 0000000..d27e6de --- /dev/null +++ b/Trivac/src/FLDTHR.f @@ -0,0 +1,300 @@ +*DECK FLDTHR + SUBROUTINE FLDTHR(IPTRK,IPSYS,IPFLUX,LADJ,LL4,ITY,NUN,NGRP,ICL1, + 1 ICL2,IMPX,NADI,NSTARD,MAXINR,EPSINR,ITER,TKT,TKB,GRAD1) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Perform thermal (up-scattering) iterations in Trivac. +* +*Copyright: +* Copyright (C) 2002 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. +* IPFLUX L_FLUX pointer to the solution. +* LADJ flag set to .TRUE. for adjoint solution acceleration. +* 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 iterations in one cycle of the inverse power +* method. +* ICL2 number of accelerated iterations in one cycle. +* IMPX print parameter (set to 0 for no printing). +* NADI number of inner ADI iterations per outer iteration. +* NSTARD number of restarting iterations with GMRES. +* MAXINR maximum number of thermal iterations. +* EPSINR thermal iteration epsilon. +* +*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. +* GRAD1 delta flux for this outer iteration. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPTRK,IPSYS,IPFLUX + INTEGER LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,NADI,NSTARD,MAXINR,ITER + REAL EPSINR,TKT,TKB,GRAD1(NUN,NGRP) + LOGICAL LADJ +*---- +* LOCAL VARIABLES +*---- + PARAMETER(NSTATE=40) + INTEGER ISTATE(NSTATE) + REAL(KIND=8) DERTOL + CHARACTER TEXT12*12,TEXT3*3 + INTERFACE + FUNCTION FLDONE_TEMPLATE(X,B,N,IPTRK,IPSYS,IPFLUX) RESULT(Y) + USE GANLIB + INTEGER, INTENT(IN) :: N + REAL(KIND=8), DIMENSION(N), INTENT(IN) :: X, B + REAL(KIND=8), DIMENSION(N) :: Y + TYPE(C_PTR) IPTRK,IPSYS,IPFLUX + END FUNCTION FLDONE_TEMPLATE + END INTERFACE + PROCEDURE(FLDONE_TEMPLATE) :: FLDONE +*---- +* ALLOCATABLE ARRAYS +*---- + REAL, DIMENSION(:), ALLOCATABLE :: W + REAL, DIMENSION(:,:), ALLOCATABLE :: GAR2 + REAL, DIMENSION(:,:,:), ALLOCATABLE :: WORK + REAL, DIMENSION(:), POINTER :: AGAR + REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: DWORK1,DWORK2 + TYPE(C_PTR) AGAR_PTR +*---- +* SCRATCH STORAGE ALLOCATION +*---- + IF(MAXINR.EQ.0) RETURN + ALLOCATE(GAR2(NUN,NGRP),WORK(LL4,NGRP,3)) +* + IF(NSTARD.GT.0) CALL LCMGET(IPFLUX,'STATE-VECTOR',ISTATE) + NCTOT=ICL1+ICL2 + IF(ICL2.EQ.0) THEN + NCPTM=NCTOT+1 + ELSE + NCPTM=ICL1 + ENDIF + DO 11 IGR=1,NGRP + DO 10 I=1,LL4 + WORK(I,IGR,1)=0.0 + WORK(I,IGR,2)=0.0 + WORK(I,IGR,3)=GRAD1(I,IGR) + 10 CONTINUE + 11 CONTINUE + IGDEB=1 +*---- +* PERFORM THERMAL (UP-SCATTERING) ITERATIONS +*---- + TEXT3='NO ' + ITER=2 + DO + CALL KDRCPU(TK1) + IF(LADJ) THEN +* ADJOINT SOLUTION + DO 31 IGR=IGDEB,NGRP + WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,IGR,3), + 1 GAR2(1,IGR)) + DO 30 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(W(LL4)) + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,JGR,3), + 1 W(1)) + DO 15 I=1,LL4 + GAR2(I,IGR)=GAR2(I,IGR)-W(I) + 15 CONTINUE + DEALLOCATE(W) + ELSE + CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR) + CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /)) + DO 20 I=1,ILONG + GAR2(I,IGR)=GAR2(I,IGR)-AGAR(I)*WORK(I,JGR,3) + 20 CONTINUE + ENDIF + 30 CONTINUE + 31 CONTINUE + DO 61 IGR=NGRP,IGDEB,-1 + DO 50 JGR=NGRP,IGR+1,-1 + WRITE(TEXT12,'(1HA,2I3.3)') JGR,IGR + CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM) + IF(ILONG.EQ.0) GO TO 50 + IF(ITY.EQ.13) THEN + ALLOCATE(W(LL4)) + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,JGR),W(1)) + DO 35 I=1,LL4 + GAR2(I,IGR)=GAR2(I,IGR)+W(I) + 35 CONTINUE + DEALLOCATE(W) + ELSE + CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR) + CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /)) + DO 40 I=1,ILONG + GAR2(I,IGR)=GAR2(I,IGR)+AGAR(I)*GAR2(I,JGR) + 40 CONTINUE + ENDIF + 50 CONTINUE + CALL KDRCPU(TK2) + TKB=TKB+(TK2-TK1) +* + CALL KDRCPU(TK1) + IF(NSTARD.EQ.0) THEN + WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR + CALL FLDADI(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,IGR),NADI) + JTER=NADI + ELSE +* use a GMRES solution of the linear system + DERTOL=EPSINR + ISTATE(39)=IGR + CALL LCMPUT(IPFLUX,'STATE-VECTOR',NSTATE,1,ISTATE) + ALLOCATE(DWORK1(LL4),DWORK2(LL4)) + DWORK1(:LL4)=GAR2(:LL4,IGR) ! source + DWORK2(:LL4)=WORK(:LL4,IGR,3) ! estimate of the flux + CALL FLDMRA(DWORK1,FLDONE,LL4,DERTOL,NSTARD,NADI,IMPX, + 1 IPTRK,IPSYS,IPFLUX,DWORK2,JTER) + GAR2(:LL4,IGR)=REAL(DWORK2(:LL4)) + DEALLOCATE(DWORK2,DWORK1) + ENDIF + DO 60 I=1,LL4 + WORK(I,IGR,1)=WORK(I,IGR,2) + WORK(I,IGR,2)=WORK(I,IGR,3) + WORK(I,IGR,3)=GRAD1(I,IGR)+(WORK(I,IGR,2)-GAR2(I,IGR)) + 60 CONTINUE + 61 CONTINUE + ELSE +* DIRECT SOLUTION + DO 81 IGR=IGDEB,NGRP + WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,IGR,3), + 1 GAR2(1,IGR)) + DO 80 JGR=1,NGRP + IF(JGR.EQ.IGR) GO TO 80 + WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR + CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM) + IF(ILONG.EQ.0) GO TO 80 + IF(ITY.EQ.13) THEN + ALLOCATE(W(LL4)) + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,JGR,3), + 1 W(1)) + DO 65 I=1,LL4 + GAR2(I,IGR)=GAR2(I,IGR)-W(I) + 65 CONTINUE + DEALLOCATE(W) + ELSE + CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR) + CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /)) + DO 70 I=1,ILONG + GAR2(I,IGR)=GAR2(I,IGR)-AGAR(I)*WORK(I,JGR,3) + 70 CONTINUE + ENDIF + 80 CONTINUE + 81 CONTINUE + DO 115 IGR=IGDEB,NGRP + DO 100 JGR=1,IGR-1 + WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR + CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM) + IF(ILONG.EQ.0) GO TO 100 + IF(ITY.EQ.13) THEN + ALLOCATE(W(LL4)) + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,JGR),W(1)) + DO 85 I=1,LL4 + GAR2(I,IGR)=GAR2(I,IGR)+W(I) + 85 CONTINUE + DEALLOCATE(W) + ELSE + CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR) + CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /)) + DO 90 I=1,ILONG + GAR2(I,IGR)=GAR2(I,IGR)+AGAR(I)*GAR2(I,JGR) + 90 CONTINUE + ENDIF + 100 CONTINUE + CALL KDRCPU(TK2) + TKB=TKB+(TK2-TK1) +* + CALL KDRCPU(TK1) + WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR + IF(NSTARD.EQ.0) THEN + WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR + CALL FLDADI(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,IGR),NADI) + JTER=NADI + ELSE +* use a GMRES solution of the linear system + DERTOL=EPSINR + ISTATE(39)=IGR + CALL LCMPUT(IPFLUX,'STATE-VECTOR',NSTATE,1,ISTATE) + ALLOCATE(DWORK1(LL4),DWORK2(LL4)) + DWORK1(:LL4)=GAR2(:LL4,IGR) ! source + DWORK2(:LL4)=WORK(:LL4,IGR,3) ! estimate of the flux + CALL FLDMRA(DWORK1,FLDONE,LL4,DERTOL,NSTARD,NADI,IMPX, + 1 IPTRK,IPSYS,IPFLUX,DWORK2,JTER) + GAR2(:LL4,IGR)=REAL(DWORK2(:LL4)) + DEALLOCATE(DWORK2,DWORK1) + ENDIF + DO 110 I=1,LL4 + WORK(I,IGR,1)=WORK(I,IGR,2) + WORK(I,IGR,2)=WORK(I,IGR,3) + WORK(I,IGR,3)=GRAD1(I,IGR)+(WORK(I,IGR,2)-GAR2(I,IGR)) + 110 CONTINUE + 115 CONTINUE + ENDIF + IF(MOD(ITER-2,NCTOT).GE.NCPTM) THEN + CALL FLD2AC(NGRP,LL4,IGDEB,WORK,ZMU) + ELSE + ZMU=1.0 + ENDIF + IGDEBO=IGDEB + DO 130 IGR=IGDEBO,NGRP + GINN=0.0 + FINN=0.0 + DO 120 I=1,LL4 + GINN=MAX(GINN,ABS(WORK(I,IGR,2)-WORK(I,IGR,3))) + FINN=MAX(FINN,ABS(WORK(I,IGR,3))) + 120 CONTINUE + GINN=GINN/FINN + IF((GINN.LT.EPSINR).AND.(IGDEB.EQ.IGR)) IGDEB=IGDEB+1 + 130 CONTINUE + CALL KDRCPU(TK2) + TKT=TKT+(TK2-TK1) + IF(GINN.LT.EPSINR) TEXT3='YES' + IF(IMPX.GT.2) WRITE(6,1000) ITER,GINN,EPSINR,IGDEB,ZMU,TEXT3, + 1 JTER + IF((GINN.LT.EPSINR).OR.(ITER.EQ.MAXINR)) EXIT + ITER=ITER+1 + ENDDO +*---- +* END OF THERMAL ITERATIONS +*---- + DO 175 I=1,LL4 + DO 170 IGR=1,NGRP + GRAD1(I,IGR)=WORK(I,IGR,3) + 170 CONTINUE + 175 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(GAR2,WORK) + RETURN +* + 1000 FORMAT (10X,3HIN(,I3,6H) FLX:,5H PRC=,1P,E9.2,5H TAR=,E9.2, + 1 7H IGDEB=, I13,6H ACCE=,0P,F12.5,12H CONVERGED=,A3,6H JTER=, + 2 I4) + END -- cgit v1.2.3