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/VALU5C.f | 133 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100755 Trivac/src/VALU5C.f (limited to 'Trivac/src/VALU5C.f') diff --git a/Trivac/src/VALU5C.f b/Trivac/src/VALU5C.f new file mode 100755 index 0000000..53539a8 --- /dev/null +++ b/Trivac/src/VALU5C.f @@ -0,0 +1,133 @@ +*DECK VALU5C + SUBROUTINE VALU5C (KPMAC,NX,NUN,NMIX,X,XXX,EVT,ISS,IXLG,ITRIAL, + 1 AXY) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Interpolation of the flux distribution for nodal method in 1D. +* +*Copyright: +* Copyright (C) 2021 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 +* KPMAC group directory in the macrolib. +* NX number of elements along the X axis. +* NUN dimension of unknown array EVT. +* NMIX number of mixtures. +* X Cartesian coordinates along the X axis where the flux is +* interpolated. +* XXX Cartesian coordinates along the X axis. +* EVT reconstruction coefficients of the flux. +* ISS mixture index assigned to each element. +* IXLG number of interpolated points according to X. +* ITRIAL type of expansion functions in the nodal calculation +* (=0: CMFD; =1: polynomial NEM; =2: hyperbolic NEM). +* +*Parameters: output +* AXY interpolated fluxes. +* +*---------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) KPMAC + INTEGER NX,NUN,NMIX,ISS(NX),IXLG,ITRIAL + REAL X(IXLG),XXX(NX+1),EVT(NUN),AXY(IXLG) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION WORK1(4,5),WORK2(2,3) + DOUBLE PRECISION GAR,ETA,ALP1,COEF,U + REAL, ALLOCATABLE, DIMENSION(:) :: DIFF,SIGR,SIGW +*---- +* RECOVER REMOVAL CROSS SECTIONS AND DIFFUSION COEFFICIENTS +*---- + ALLOCATE(DIFF(NMIX),SIGR(NMIX),SIGW(NMIX)) + CALL LCMGET(KPMAC,'NTOT0',SIGR) + CALL LCMGET(KPMAC,'SIGW00',SIGW) + CALL LCMGET(KPMAC,'DIFF',DIFF) + SIGR(:)=SIGR(:)-SIGW(:) +*---- +* PERFORM INTERPOLATION +*---- + DO I=1,IXLG + ABSC=X(I) + GAR=0.0D0 +* +* Find the node index containing the interpolation point + IS=0 + DO KEL=1,NX + IS=KEL + IF((ABSC.GE.XXX(KEL)).AND.(ABSC.LE.XXX(KEL+1))) GO TO 10 + ENDDO + CALL XABORT('VALU5C: WRONG INTERPOLATION.') + 10 IBM=ISS(IS) + IF(IBM.EQ.0) GO TO 100 + ETA=(XXX(IS+1)-XXX(IS))*SQRT(SIGR(IBM)/DIFF(IBM)) + ALP1=ETA*COSH(ETA/2.0)-2.0*SINH(ETA/2.0) + COEF=DIFF(IBM)/(XXX(IS+1)-XXX(IS)) + U=(ABSC-XXX(IS))/(XXX(IS+1)-XXX(IS))-0.5 + IF(ITRIAL.EQ.0) THEN + WORK2(1,1)=COEF + WORK2(1,2)=-3.0*COEF + WORK2(1,3)=EVT(3*NX+IS) + WORK2(2,1)=COEF + WORK2(2,2)=3.0*COEF + WORK2(2,3)=EVT(3*NX+IS+1) + CALL ALSBD(3,1,WORK2(1,1),IER,3) + IF(IER.NE.0) CALL XABORT('VALU5C: SINGULAR MATRIX(1).') + GAR=EVT(IS)+WORK2(1,3)*U+WORK2(2,3)*(3.0*U**2-1.0/4.0) + ELSE + WORK1(:,:)=0.0 + WORK1(1,1)=-0.5 + WORK1(1,2)=0.5 + WORK1(1,5)=EVT(NX+IS)-EVT(IS) + WORK1(2,1)=0.5 + WORK1(2,2)=0.5 + WORK1(2,5)=EVT(2*NX+IS)-EVT(IS) + WORK1(3,1)=-COEF + WORK1(3,2)=3.0*COEF + WORK1(3,5)=EVT(3*NX+IS) + WORK1(4,1)=-COEF + WORK1(4,2)=-3.0*COEF + WORK1(4,5)=EVT(3*NX+IS+1) + IF(ITRIAL.EQ.1) THEN + WORK1(3,3)=-0.5*COEF + WORK1(3,4)=0.2*COEF + WORK1(4,3)=-0.5*COEF + WORK1(4,4)=-0.2*COEF + ELSE + WORK1(1,3)=-SINH(ETA/2.0) + WORK1(1,4)=ALP1/ETA + WORK1(2,3)=SINH(ETA/2.0) + WORK1(2,4)=ALP1/ETA + WORK1(3,3)=-COEF*ETA*COSH(ETA/2.0) + WORK1(3,4)=COEF*ETA*SINH(ETA/2.0) + WORK1(4,3)=-COEF*ETA*COSH(ETA/2.0) + WORK1(4,4)=-COEF*ETA*SINH(ETA/2.0) + ENDIF + CALL ALSBD(4,1,WORK1(1,1),IER,4) + IF(IER.NE.0) CALL XABORT('VALU5C: SINGULAR MATRIX(2).') + GAR=EVT(IS)+WORK1(1,5)*U+WORK1(2,5)*(3.0*U**2-1.0/4.0) + IF(ITRIAL.EQ.1) THEN + GAR=GAR+WORK1(3,5)*(U**2-0.25)*U+ + 1 WORK1(4,5)*(U**2-0.25)*(U**2-0.05) + ELSE + GAR=GAR+WORK1(3,5)*SINH(ETA*U)+ + 1 WORK1(4,5)*(COSH(ETA*U)-2.0*SINH(ETA/2.0)/ETA) + ENDIF + ENDIF + 100 AXY(I)=REAL(GAR) + ENDDO + DEALLOCATE(SIGW,SIGR,DIFF) + RETURN + END -- cgit v1.2.3