summaryrefslogtreecommitdiff
path: root/Trivac/src/VALU5B.f
diff options
context:
space:
mode:
Diffstat (limited to 'Trivac/src/VALU5B.f')
-rwxr-xr-xTrivac/src/VALU5B.f342
1 files changed, 342 insertions, 0 deletions
diff --git a/Trivac/src/VALU5B.f b/Trivac/src/VALU5B.f
new file mode 100755
index 0000000..fe61753
--- /dev/null
+++ b/Trivac/src/VALU5B.f
@@ -0,0 +1,342 @@
+*DECK VALU5B
+ SUBROUTINE VALU5B (KPMAC,NX,NY,LL4F,LL4X,NUN,NMIX,X,Y,XXX,YYY,
+ 1 EVT,ISS,KFLX,KN,IXLG,IYLG,ICORN,AXY)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Interpolation of the flux distribution for nodal method in 2D.
+*
+*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.
+* NY number of elements along the Y axis.
+* LL4F number of averaged flux unknowns.
+* LL4X number of X-directed net currents.
+* NUN dimension of unknown array EVT.
+* NMIX number of mixtures.
+* X Cartesian coordinates along the X axis where the flux is
+* interpolated.
+* Y Cartesian coordinates along the Y axis where the flux is
+* interpolated.
+* XXX Cartesian coordinates along the X axis.
+* YYY Cartesian coordinates along the Y axis.
+* EVT reconstruction coefficients of the flux.
+* ISS mixture index assigned to each element.
+* KFLX correspondence between local and global numbering.
+* KN element-ordered interface net current unknown list.
+* IXLG number of interpolated points according to X.
+* IYLG number of interpolated points according to Y.
+* ICORN flag to activate corner flux correction (0/1: OFF/ON).
+*
+*Parameters: output
+* AXY interpolated fluxes.
+*
+*----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) KPMAC
+ INTEGER NX,NY,LL4F,LL4X,NUN,NMIX,ISS(NX*NY),KFLX(NX*NY),
+ 1 KN(6,NX,NY),IXLG,IYLG,ICORN
+ REAL X(IXLG),Y(IYLG),XXX(NX+1),YYY(NY+1),EVT(NUN),AXY(IXLG,IYLG)
+*----
+* LOCAL VARIABLES
+*----
+ DOUBLE PRECISION WORK1(4,5),FC2(4)
+ DOUBLE PRECISION GAR,COEFX,COEFY,U,V,P2U,P2V
+ LOGICAL LOGC1,LOGC2
+*----
+* ALLOCATABLE ARRAYS
+*----
+ REAL, ALLOCATABLE, DIMENSION(:) :: DIFF
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: FCORN,DELC
+*----
+* RECOVER DIFFUSION COEFFICIENTS
+*----
+ ALLOCATE(DIFF(NMIX))
+ CALL LCMGET(KPMAC,'DIFF',DIFF)
+*----
+* COMPUTE CORNER FLUXES
+*----
+ ALLOCATE(DELC(4,NX,NY))
+ DELC(:4,:NX,:NY)=0.D0
+ IF(ICORN==1) THEN
+ ALLOCATE(FCORN(4,NX,NY))
+ FCORN(:4,:NX,:NY)=0.D0
+ DO JS=1,NY
+ DO IS=1,NX
+ IEL=(JS-1)*NX+IS
+ IND1=KFLX(IEL)
+ IF(IND1.EQ.0) CYCLE
+ IBM=ISS(IEL)
+ IF(IBM.LE.0) CYCLE
+ JXM=KN(1,IS,JS) ; JXP=KN(2,IS,JS)
+ JYM=KN(3,IS,JS) ; JYP=KN(4,IS,JS)
+ COEFX=DIFF(IBM)/(XXX(IS+1)-XXX(IS))
+ COEFY=DIFF(IBM)/(YYY(JS+1)-YYY(JS))
+*
+ WORK1(:,:)=0.0
+ WORK1(1,1)=-0.5
+ WORK1(1,2)=0.5
+ WORK1(1,5)=EVT(LL4F+IND1)-EVT(IND1)
+ WORK1(2,1)=0.5
+ WORK1(2,2)=0.5
+ WORK1(2,5)=EVT(2*LL4F+IND1)-EVT(IND1)
+ WORK1(3,1)=-COEFX
+ WORK1(3,2)=3.0*COEFX
+ IF(JXM.NE.0) WORK1(3,5)=EVT(5*LL4F+JXM)
+ WORK1(4,1)=-COEFX
+ WORK1(4,2)=-3.0*COEFX
+ IF(JXP.NE.0) WORK1(4,5)=EVT(5*LL4F+JXP)
+ WORK1(3,3)=-0.5*COEFX
+ WORK1(3,4)=0.2*COEFX
+ WORK1(4,3)=-0.5*COEFX
+ WORK1(4,4)=-0.2*COEFX
+ CALL ALSBD(4,1,WORK1,IER,4)
+ IF(IER.NE.0) CALL XABORT('VALU5B: SINGULAR MATRIX(1).')
+ DO IC=1,4
+ SELECT CASE(IC)
+ CASE(1,3)
+ U=-0.5
+ CASE DEFAULT
+ U=0.5
+ END SELECT
+ GAR=EVT(IND1)+WORK1(1,5)*U+WORK1(2,5)*(3.0*U**2-0.25)
+ GAR=GAR+WORK1(3,5)*(U**2-0.25)*U+WORK1(4,5)*(U**2-0.25)*
+ 1 (U**2-0.05)
+ FCORN(IC,IS,JS)=GAR
+ ENDDO
+*
+ WORK1(:,:)=0.0
+ WORK1(1,1)=-0.5
+ WORK1(1,2)=0.5
+ WORK1(1,5)=EVT(3*LL4F+IND1)-EVT(IND1)
+ WORK1(2,1)=0.5
+ WORK1(2,2)=0.5
+ WORK1(2,5)=EVT(4*LL4F+IND1)-EVT(IND1)
+ WORK1(3,1)=-COEFY
+ WORK1(3,2)=3.0*COEFY
+ IF(JYM.NE.0) WORK1(3,5)=EVT(5*LL4F+LL4X+JYM)
+ WORK1(4,1)=-COEFY
+ WORK1(4,2)=-3.0*COEFY
+ IF(JYP.NE.0) WORK1(4,5)=EVT(5*LL4F+LL4X+JYP)
+ WORK1(3,3)=-0.5*COEFY
+ WORK1(3,4)=0.2*COEFY
+ WORK1(4,3)=-0.5*COEFY
+ WORK1(4,4)=-0.2*COEFY
+ CALL ALSBD(4,1,WORK1,IER,4)
+ IF(IER.NE.0) CALL XABORT('VALU5B: SINGULAR MATRIX(2).')
+ DO IC=1,4
+ SELECT CASE(IC)
+ CASE(1,2)
+ V=-0.5
+ CASE DEFAULT
+ V=0.5
+ END SELECT
+ GAR=FCORN(IC,IS,JS)+WORK1(1,5)*V+WORK1(2,5)*
+ 1 (3.0*V**2-0.25)
+ GAR=GAR+WORK1(3,5)*(V**2-0.25)*V+WORK1(4,5)*(V**2-0.25)*
+ 1 (V**2-0.05)
+ FCORN(IC,IS,JS)=GAR
+ ENDDO
+ ENDDO
+ ENDDO
+ DO JS=1,NY
+ DO IS=1,NX
+ IEL=(JS-1)*NX+IS
+ IND1=KFLX(IEL)
+ IF(IND1.EQ.0) CYCLE
+ ! corner 1
+ NB=1; GAR=FCORN(1,IS,JS)
+ LOGC1=(IS>1) ; LOGC2=(JS>1)
+ IF(LOGC2) LOGC2=(KFLX((JS-2)*NX+IS)>0)
+ IF(LOGC1) THEN
+ IF(KFLX((JS-1)*NX+IS-1)>0) THEN
+ NB=NB+1 ;GAR=GAR+FCORN(2,IS-1,JS)
+ ENDIF
+ ENDIF
+ IF(LOGC2) THEN
+ IF(KFLX((JS-2)*NX+IS)>0) THEN
+ NB=NB+1 ;GAR=GAR+FCORN(3,IS,JS-1)
+ ENDIF
+ ENDIF
+ IF(LOGC1.AND.LOGC2) THEN
+ IF(KFLX((JS-2)*NX+IS-1)>0) THEN
+ NB=NB+1 ;GAR=GAR+FCORN(4,IS-1,JS-1)
+ ENDIF
+ ENDIF
+ FC2(1)=GAR/REAL(NB)-FCORN(1,IS,JS)
+ ! corner 2
+ NB=1 ;GAR=FCORN(2,IS,JS)
+ LOGC1=(IS<NX) ; LOGC2=(JS>1)
+ IF(LOGC1) THEN
+ IF(KFLX((JS-1)*NX+IS+1)>0) THEN
+ NB=NB+1 ;GAR=GAR+FCORN(1,IS+1,JS)
+ ENDIF
+ ENDIF
+ IF(LOGC2) THEN
+ IF(KFLX((JS-2)*NX+IS)>0) THEN
+ NB=NB+1 ;GAR=GAR+FCORN(4,IS,JS-1)
+ ENDIF
+ ENDIF
+ IF(LOGC1.AND.LOGC2) THEN
+ IF(KFLX((JS-2)*NX+IS+1)>0) THEN
+ NB=NB+1 ;GAR=GAR+FCORN(3,IS+1,JS-1)
+ ENDIF
+ ENDIF
+ FC2(2)=GAR/REAL(NB)-FCORN(2,IS,JS)
+ ! corner 3
+ NB=1 ; GAR=FCORN(3,IS,JS)
+ LOGC1=(IS>1) ; LOGC2=(JS<NY)
+ IF(LOGC1) THEN
+ IF(KFLX((JS-1)*NX+IS-1)>0) THEN
+ NB=NB+1 ; GAR=GAR+FCORN(4,IS-1,JS)
+ ENDIF
+ ENDIF
+ IF(LOGC2) THEN
+ IF(KFLX(JS*NX+IS)>0) THEN
+ NB=NB+1 ; GAR=GAR+FCORN(1,IS,JS+1)
+ ENDIF
+ ENDIF
+ IF(LOGC1.AND.LOGC2) THEN
+ IF(KFLX(JS*NX+IS-1)>0) THEN
+ NB=NB+1 ; GAR=GAR+FCORN(2,IS-1,JS+1)
+ ENDIF
+ ENDIF
+ FC2(3)=GAR/REAL(NB)-FCORN(3,IS,JS)
+ ! corner 4
+ NB=1
+ GAR=FCORN(4,IS,JS)
+ LOGC1=(IS<NX)
+ IF(LOGC1) LOGC1=(KFLX((JS-1)*NX+IS+1)>0)
+ LOGC2=(JS<NY)
+ IF(LOGC2) LOGC2=(KFLX(JS*NX+IS)>0)
+ IF(LOGC1) THEN
+ IF(KFLX((JS-1)*NX+IS+1)>0) THEN
+ NB=NB+1 ; GAR=GAR+FCORN(3,IS+1,JS)
+ ENDIF
+ ENDIF
+ IF(LOGC2) THEN
+ IF(KFLX(JS*NX+IS)>0) THEN
+ NB=NB+1 ; GAR=GAR+FCORN(2,IS,JS+1)
+ ENDIF
+ ENDIF
+ IF(LOGC1.AND.LOGC2) THEN
+ IF(KFLX(JS*NX+IS+1)>0) THEN
+ NB=NB+1 ; GAR=GAR+FCORN(1,IS+1,JS+1)
+ ENDIF
+ ENDIF
+ FC2(4)=GAR/REAL(NB)-FCORN(4,IS,JS)
+ ! polynomial coefficients of correction terms
+ DELC(1,IS,JS)= FC2(1)-FC2(2)-FC2(3)+FC2(4)
+ DELC(2,IS,JS)=-FC2(1)-FC2(2)+FC2(3)+FC2(4)
+ DELC(3,IS,JS)=-FC2(1)+FC2(2)-FC2(3)+FC2(4)
+ DELC(4,IS,JS)= FC2(1)+FC2(2)+FC2(3)+FC2(4)
+ ENDDO
+ ENDDO
+ DEALLOCATE(FCORN)
+ ENDIF
+*----
+* PERFORM INTERPOLATION
+*----
+ DO J=1,IYLG
+ ORDO=Y(J)
+ DO I=1,IXLG
+ ABSC=X(I)
+ GAR=0.0D0
+*
+* Find the node index containing the interpolation point
+ IS=0; JS=0
+ DO L=1,NX
+ IS=L
+ IF((ABSC.GE.XXX(L)).AND.(ABSC.LE.XXX(L+1))) GO TO 10
+ ENDDO
+ CALL XABORT('VALU5B: WRONG INTERPOLATION(1).')
+ 10 DO L=1,NY
+ JS=L
+ IF((ORDO.GE.YYY(L)).AND.(ORDO.LE.YYY(L+1))) GO TO 20
+ ENDDO
+ CALL XABORT('VALU5B: WRONG INTERPOLATION(2).')
+ 20 IEL=(JS-1)*NX+IS
+ IND1=KFLX(IEL)
+ IF(IND1.EQ.0) GO TO 30
+ IBM=ISS(IEL)
+ IF(IBM.LE.0) GO TO 30
+ JXM=KN(1,IS,JS) ; JXP=KN(2,IS,JS)
+ JYM=KN(3,IS,JS) ; JYP=KN(4,IS,JS)
+ COEFX=DIFF(IBM)/(XXX(IS+1)-XXX(IS))
+ COEFY=DIFF(IBM)/(YYY(JS+1)-YYY(JS))
+ U=(ABSC-XXX(IS))/(XXX(IS+1)-XXX(IS))-0.5
+ V=(ORDO-YYY(JS))/(YYY(JS+1)-YYY(JS))-0.5
+ GAR=EVT(IND1)
+*
+ WORK1(:,:)=0.0
+ WORK1(1,1)=-0.5
+ WORK1(1,2)=0.5
+ WORK1(1,5)=EVT(LL4F+IND1)-EVT(IND1)
+ WORK1(2,1)=0.5
+ WORK1(2,2)=0.5
+ WORK1(2,5)=EVT(2*LL4F+IND1)-EVT(IND1)
+ WORK1(3,1)=-COEFX
+ WORK1(3,2)=3.0*COEFX
+ IF(JXM.NE.0) WORK1(3,5)=EVT(5*LL4F+JXM)
+ WORK1(4,1)=-COEFX
+ WORK1(4,2)=-3.0*COEFX
+ IF(JXP.NE.0) WORK1(4,5)=EVT(5*LL4F+JXP)
+ WORK1(3,3)=-0.5*COEFX
+ WORK1(3,4)=0.2*COEFX
+ WORK1(4,3)=-0.5*COEFX
+ WORK1(4,4)=-0.2*COEFX
+ CALL ALSBD(4,1,WORK1,IER,4)
+ IF(IER.NE.0) CALL XABORT('VALU5B: SINGULAR MATRIX(3).')
+ GAR=GAR+WORK1(1,5)*U+WORK1(2,5)*(3.0*U**2-0.25)
+ GAR=GAR+WORK1(3,5)*(U**2-0.25)*U+WORK1(4,5)*(U**2-0.25)*
+ 1 (U**2-0.05)
+*
+ WORK1(:,:)=0.0
+ WORK1(1,1)=-0.5
+ WORK1(1,2)=0.5
+ WORK1(1,5)=EVT(3*LL4F+IND1)-EVT(IND1)
+ WORK1(2,1)=0.5
+ WORK1(2,2)=0.5
+ WORK1(2,5)=EVT(4*LL4F+IND1)-EVT(IND1)
+ WORK1(3,1)=-COEFY
+ WORK1(3,2)=3.0*COEFY
+ IF(JYM.NE.0) WORK1(3,5)=EVT(5*LL4F+LL4X+JYM)
+ WORK1(4,1)=-COEFY
+ WORK1(4,2)=-3.0*COEFY
+ IF(JYP.NE.0) WORK1(4,5)=EVT(5*LL4F+LL4X+JYP)
+ WORK1(3,3)=-0.5*COEFY
+ WORK1(3,4)=0.2*COEFY
+ WORK1(4,3)=-0.5*COEFY
+ WORK1(4,4)=-0.2*COEFY
+ CALL ALSBD(4,1,WORK1,IER,4)
+ IF(IER.NE.0) CALL XABORT('VALU5B: SINGULAR MATRIX(4).')
+ GAR=GAR+WORK1(1,5)*V+WORK1(2,5)*(3.0*V**2-0.25)
+ GAR=GAR+WORK1(3,5)*(V**2-0.25)*V+WORK1(4,5)*(V**2-0.25)*
+ 1 (V**2-0.05)
+*
+ IF(ICORN==1) THEN
+ ! perform interpolation of corner flux correction
+ P2U=3.0*U**2-0.25 ; P2V=3.0*V**2-0.25
+ GAR=GAR+DELC(1,IS,JS)*U*V + DELC(2,IS,JS)*P2U*V+
+ 1 DELC(3,IS,JS)*U*P2V + DELC(4,IS,JS)*P2U*P2V
+ ENDIF
+ 30 AXY(I,J)=REAL(GAR)
+ ENDDO
+ ENDDO
+ DEALLOCATE(DELC,DIFF)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:26:40 +0000