summaryrefslogtreecommitdiff
path: root/Trivac/src/VALU5C.f
blob: 53539a8b598d1e8a9270f5edf30f7268ce6b3355 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
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