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 --- Dragon/src/MCGCDD.f | 176 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 Dragon/src/MCGCDD.f (limited to 'Dragon/src/MCGCDD.f') diff --git a/Dragon/src/MCGCDD.f b/Dragon/src/MCGCDD.f new file mode 100644 index 0000000..97c2158 --- /dev/null +++ b/Dragon/src/MCGCDD.f @@ -0,0 +1,176 @@ +*DECK MCGCDD + SUBROUTINE MCGCDD(IPRINT,IPMACR,II,NG,NGEFF,NGIND,NCONV,M,NLONG, + 1 NUN,NREG,LC,LFORW,PACA,NZON,KEYFLX,KEYCUR,IPERM, + 2 IM,MCU,JU,EPSINT,MAXINT,FLUX,QFR,DIAGQ,CQ,DIAGF, + 3 CF,ILUDF,ILUCF,LC0,IM0,MCU0) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solution of the CDD equations (ACA method) for a synthetic diffusion +* flux calculation. +* +*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): I. Suslov and R. Le Tellier +* +*Parameters: input +* IPRINT print parameter. +* IPMACR pointer to the macrolib LCM object ('GROUP' directory) +* II group being processed (II <= NGEFF). +* NG number of groups. +* NGEFF number of groups to process. +* NGIND index of the groups to process. +* NCONV logical array of convergence status for each group (.TRUE. +* not converged). +* M number of material mixtures. +* NLONG size of the corrective system. +* NUN number of unknowns per group. +* NREG number of volume regions. +* LC dimension of profiled matrices MCU and CQ. +* LFORW flag set to .false. to transpose the coefficient matrix. +* PACA type of preconditioner to solve the ACA corrective system. +* NZON index-number of the mixture type assigned to each volume. +* KEYFLX position of flux elements in FLUX vector. +* KEYCUR position of current elements in FLUX vector. +* IPERM permutation array. +* IM used in cdd acceleration. +* MCU used in cdd acceleration. +* JU used in ACA acceleration for ilu0. +* EPSINT stopping criterion for BICGSTAB in ACA resolution. +* MAXINT maximum number of iterations allowed for BICGSTAB in ACA +* resolution. +* QFR source vector. +* DIAGQ used in cdd acceleration. +* CQ used in cdd acceleration. +* DIAGF used in cdd acceleration. +* CF used in cdd acceleration. +* ILUDF used in cdd acceleration. +* ILUCF used in cdd acceleration. +* LC0 used in ILU0-ACA acceleration. +* IM0 used in ILU0-ACA acceleration. +* MCU0 used in ILU0-ACA acceleration. +* +*Parameters: output +* FLUX zonal scalar flux. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPMACR + INTEGER IPRINT,II,NG,NGEFF,NGIND(NGEFF),M,NLONG,NUN,NREG,LC,PACA, + 1 NZON(NLONG),KEYFLX(NREG),KEYCUR(*),IPERM(NLONG),IM(NLONG+1), + 2 MCU(LC),JU(NLONG),MAXINT,LC0,IM0(*),MCU0(*) + REAL EPSINT,FLUX(NUN),QFR(NUN,NGEFF),DIAGQ(NLONG),CQ(LC), + 1 DIAGF(NLONG),CF(LC),ILUDF(NLONG),ILUCF(LC) + LOGICAL LFORW,NCONV(NGEFF) +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) JPMACR,KPMACR + DOUBLE PRECISION FF +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: NJJ,IJJ,IPOS + REAL, ALLOCATABLE, DIMENSION(:) :: XSCAT + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: AR,PHI +*---- +* SCRATCH STORAGE ALLOCATION +*---- + IF(C_ASSOCIATED(IPMACR)) THEN + JPMACR=LCMGID(IPMACR,'GROUP') + ALLOCATE(NJJ(0:M),IJJ(0:M),IPOS(0:M),XSCAT(0:M*NG)) + ELSE + JPMACR=C_NULL_PTR + ENDIF + ALLOCATE(PHI(NLONG),AR(NLONG)) +*---- +* CONSTRUCT RHS OF THE CDD SYSTEM +*---- + DO I=1,NLONG + J=IPERM(I) + IF(NZON(J).GE.0) THEN + FF=DIAGQ(I)*QFR(KEYFLX(J),II) + ELSE + FF=0.0 + ENDIF + DO IL=IM(I)+1,IM(I+1) + IF(MCU(IL).GT.0) THEN + L=IPERM(MCU(IL)) + IF(NZON(L).GE.0) FF=FF+CQ(IL)*QFR(KEYFLX(L),II) + ENDIF + ENDDO + PHI(I)=FF + ENDDO +*---- +* INVERSE THE SYSTEM BY THE ITERATIVE METHOD BICGSTAB +*---- +* apply preconditioner to RHS + CALL MCGPRA(LFORW,2,PACA,.TRUE.,NLONG,LC,IM,MCU,JU,DIAGF,CF, + 1 ILUDF,ILUCF,DIAGF,AR,PHI,LC0,IM0,MCU0,CF) +* + CALL MCGABG(IPRINT,LFORW,PACA,NLONG,LC,EPSINT,MAXINT,IM,MCU,JU, + 1 DIAGF,CF,ILUDF,ILUCF,AR,PHI,1.0,LC0,IM0,MCU0) +* + IF(C_ASSOCIATED(JPMACR)) THEN +*--- +* MODIFY THE CONTRIBUTION FORM THIS GROUP TO OTHER GROUP ISOTROPIC SOURCES +* (JACOBI -> GAUSS-SEIDEL) +*--- + IG=NGIND(II) + DO JJ=1,NGEFF + IF(NCONV(JJ)) THEN + JG=NGIND(JJ) + IF(JG.GT.IG) THEN + KPMACR=LCMGIL(JPMACR,JG) + CALL LCMGET(KPMACR,'NJJS00',NJJ(1)) + CALL LCMGET(KPMACR,'IJJS00',IJJ(1)) + CALL LCMGET(KPMACR,'IPOS00',IPOS(1)) + CALL LCMGET(KPMACR,'SCAT00',XSCAT(1)) + DO 10 I=1,NLONG + J=IPERM(I) + IBM=NZON(J) + IF(IBM.GT.0) THEN + IND=KEYFLX(J) + IGG=IJJ(IBM) + DO 20 JND=1,NJJ(IBM) + IF(IG.EQ.IGG) THEN + QFR(IND,JJ)=QFR(IND,JJ)+ + 1 XSCAT(IPOS(IBM)+JND-1)*(REAL(PHI(I))-FLUX(IND)) + GOTO 10 + ENDIF + IGG=IGG-1 + 20 CONTINUE + ENDIF + 10 CONTINUE + ENDIF + ENDIF + ENDDO + ENDIF +*--- +* REORDER FLUX VECTOR +*--- + DO I=1,NLONG + J=IPERM(I) + IF(NZON(J).GE.0) THEN + FLUX(KEYFLX(J))=REAL(PHI(I)) + ELSE + FLUX(KEYCUR(J-NREG))=REAL(PHI(I)) + ENDIF + ENDDO +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(AR,PHI) + IF(C_ASSOCIATED(IPMACR)) DEALLOCATE(XSCAT,IPOS,IJJ,NJJ) + RETURN + END -- cgit v1.2.3