← Back to branch summary

~fluidity-core/fluidity/exorcised

« back to all changes in this revision

Viewing changes to multiphase/src/multi_compositional.F90

  • Committer: saunde01
  • Date: 2011-03-23 13:15:31 UTC
  • Revision ID: svn-v4:5bf5533e-7014-46e3-b1bb-cce4b9d03719:trunk:3310
As discussed in the Dev meeting and in emails, this commit renames the multiphase directory to legacy_reservoir_prototype, and adds an obvious warning to the top of the README file to warn unsuspecting users of the experimental nature of the code. Iterations on the exact wording very welcome

Show diffs side-by-side

added added

removed removed

Lines of Context:
multiphase/src/multi_compositional.F90
1
 
!    Copyright (C) 2006 Imperial College London and others.
2
 
!    
3
 
!    Please see the AUTHORS file in the main source directory for a full list
4
 
!    of copyright holders.
5
 
!
6
 
!    Prof. C Pain
7
 
!    Applied Modelling and Computation Group
8
 
!    Department of Earth Science and Engineering
9
 
!    Imperial College London
10
 
!
11
 
!    amcgsoftware@imperial.ac.uk
12
 
!    
13
 
!    This library is free software; you can redistribute it and/or
14
 
!    modify it under the terms of the GNU Lesser General Public
15
 
!    License as published by the Free Software Foundation,
16
 
!    version 2.1 of the License.
17
 
!
18
 
!    This library is distributed in the hope that it will be useful,
19
 
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
20
 
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
21
 
!    Lesser General Public License for more details.
22
 
!
23
 
!    You should have received a copy of the GNU Lesser General Public
24
 
!    License along with this library; if not, write to the Free Software
25
 
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
26
 
!    USA
27
 
 
28
 
module Compositional_Terms
29
 
 
30
 
contains
31
 
 
32
 
 
33
 
  SUBROUTINE CALC_COMP_ABSORB( ICOMP, NCOMP, DT, ALPHA_BETA, &
34
 
       NPHASE, CV_NONODS, K_COMP, DENOLD, SATURAOLD, VOLFRA_PORE, COMP_ABSORB, &
35
 
       TOTELE, CV_NLOC, CV_NDGLN)
36
 
 
37
 
    ! Calculate compositional model linkage between the phase expressed in COMP_ABSORB. 
38
 
    ! Use values from the previous time step so its easier to converge. 
39
 
    ! ALPHA_BETA is the scaling coeff. of the compositional model e.g. =1.0
40
 
 
41
 
    IMPLICIT NONE
42
 
    INTEGER, intent( in ) :: ICOMP, NCOMP, NPHASE, CV_NONODS, TOTELE, CV_NLOC
43
 
    REAL, intent( in ) :: DT, ALPHA_BETA
44
 
    REAL, DIMENSION( NCOMP, NPHASE, NPHASE ), intent( in ) :: K_COMP
45
 
    REAL, DIMENSION( CV_NONODS * NPHASE ), intent( in ) :: DENOLD, SATURAOLD
46
 
    REAL, DIMENSION( TOTELE ), intent( in ) :: VOLFRA_PORE
47
 
    REAL, DIMENSION( CV_NONODS, NPHASE, NPHASE  ), intent( inout ) :: COMP_ABSORB
48
 
    INTEGER, DIMENSION( TOTELE * CV_NLOC ), intent( in ) :: CV_NDGLN
49
 
 
50
 
    ! Local Variables
51
 
    INTEGER :: IPHASE, JPHASE, ELE, CV_ILOC, CV_NOD, JCOMP
52
 
    REAL, DIMENSION( : ), allocatable :: ALPHA, SUM_NOD, VOLFRA_PORE_NOD
53
 
 
54
 
    ALLOCATE( ALPHA( CV_NONODS ))
55
 
    ALLOCATE( SUM_NOD( CV_NONODS ))
56
 
    ALLOCATE( VOLFRA_PORE_NOD( CV_NONODS ))
57
 
 
58
 
    ! Determine a node-wise representation of porosity VOLFRA_PORE_NOD.      
59
 
    SUM_NOD = 0.0
60
 
    VOLFRA_PORE_NOD = 0.0
61
 
 
62
 
    DO ELE = 1, TOTELE
63
 
       DO CV_ILOC = 1, CV_NLOC
64
 
          CV_NOD = CV_NDGLN( ( ELE - 1 ) * CV_NLOC + CV_ILOC ) 
65
 
          SUM_NOD( CV_NOD ) = SUM_NOD( CV_NOD ) + 1.0
66
 
          VOLFRA_PORE_NOD( CV_NOD ) = VOLFRA_PORE_NOD( CV_NOD ) + VOLFRA_PORE( ELE )
67
 
       END DO
68
 
    END DO
69
 
    VOLFRA_PORE_NOD = VOLFRA_PORE_NOD / SUM_NOD
70
 
 
71
 
    COMP_ABSORB = 0.0
72
 
 
73
 
    DO IPHASE = 1, NPHASE
74
 
       DO JPHASE = IPHASE + 1, NPHASE, 1
75
 
          ALPHA( 1 : CV_NONODS ) = ALPHA_BETA * VOLFRA_PORE_NOD( 1 : CV_NONODS ) * &
76
 
               ( SATURAOLD( 1 + ( IPHASE - 1 ) * CV_NONODS : IPHASE * CV_NONODS ) &
77
 
               * DENOLD( 1 +( IPHASE - 1 ) * CV_NONODS: IPHASE * CV_NONODS ) / K_COMP( ICOMP, IPHASE, JPHASE ) &
78
 
               + SATURAOLD( 1 + ( JPHASE - 1 ) * CV_NONODS : JPHASE * CV_NONODS ) * DENOLD( 1 + ( JPHASE - 1 ) * &
79
 
               CV_NONODS : JPHASE * CV_NONODS )) / DT
80
 
write(357,*) 'iphase, jphase, alpha1:',iphase, jphase, alpha 
81
 
 
82
 
          COMP_ABSORB( : , IPHASE, IPHASE ) = COMP_ABSORB( : , IPHASE, IPHASE ) &
83
 
               + ALPHA( : ) * K_COMP( ICOMP, IPHASE, JPHASE )
84
 
          COMP_ABSORB( : , IPHASE, JPHASE ) = -ALPHA( : ) 
85
 
       END DO
86
 
    END DO
87
 
 
88
 
    DO IPHASE = 1, NPHASE
89
 
       DO JPHASE = 1, IPHASE - 1
90
 
 
91
 
          ALPHA( 1 : CV_NONODS ) = ALPHA_BETA * VOLFRA_PORE_NOD( 1: CV_NONODS ) * &
92
 
               ( SATURAOLD( 1 + ( IPHASE - 1 ) * CV_NONODS : IPHASE * CV_NONODS ) * & 
93
 
               DENOLD( 1 +( IPHASE - 1 ) * CV_NONODS : IPHASE * CV_NONODS )  &
94
 
               + SATURAOLD( 1 + ( JPHASE - 1 ) * CV_NONODS : JPHASE * CV_NONODS ) * DENOLD( 1 + ( JPHASE - 1 ) * &
95
 
               CV_NONODS: JPHASE * CV_NONODS ) / K_COMP( ICOMP, JPHASE, IPHASE )) / DT
96
 
 
97
 
          COMP_ABSORB( : , IPHASE, IPHASE ) = COMP_ABSORB( : , IPHASE, IPHASE ) &
98
 
               + ALPHA( : ) 
99
 
          COMP_ABSORB( : , IPHASE, JPHASE ) = -ALPHA( : ) * K_COMP( ICOMP, JPHASE, IPHASE )
100
 
write(357,*) 'iphase, jphase, alpha2:',iphase, jphase, alpha 
101
 
 
102
 
       END DO
103
 
    END DO
104
 
 
105
 
write(357,*)'comp_ABSORB:',((comp_ABSORB(1,iphase,jphase), iphase=1,nphase),jphase=1,nphase)
106
 
write(357,*)'K_comp:'
107
 
do jcomp = 1, ncomp
108
 
do iphase = 1, nphase
109
 
write(357,*) jcomp, iphase, (K_Comp( jcomp, iphase, jphase), jphase= 1, nphase )
110
 
end do
111
 
end do
112
 
write(357,*)'K_comp:',K_comp
113
 
 
114
 
    DEALLOCATE( ALPHA )
115
 
    DEALLOCATE( SUM_NOD )
116
 
    DEALLOCATE( VOLFRA_PORE_NOD )
117
 
 
118
 
    RETURN
119
 
 
120
 
  END SUBROUTINE CALC_COMP_ABSORB
121
 
 
122
 
 
123
 
 
124
 
  SUBROUTINE CALC_COMP_DIF( NDIM, NPHASE, COMP_DIFFUSION_OPT, MAT_NONODS, &
125
 
       TOTELE, MAT_NLOC, CV_NLOC, U_NLOC, CV_SNLOC, U_SNLOC, &
126
 
       COMP_DIFFUSION, NCOMP_DIFF_COEF, COMP_DIFF_COEF, &
127
 
       X_NONODS, X, Y, Z, NU, NV, NW, U_NONODS, CV_NONODS, &
128
 
       MAT_NDGLN, U_NDGLN, X_NDGLN, &
129
 
       U_ELE_TYPE, P_ELE_TYPE ) 
130
 
    ! Calculate the diffusion coefficient COMP_DIFFUSION for current composition...
131
 
    ! based on page 136 in Reservoir-Simulation-Mathematical-Techniques-In-Oil-Recovery-(2007).pdf
132
 
    ! COMP_DIFFUSION_OPT, integer option defining diffusion coeff
133
 
    ! NCOMP_DIFF_COEF,  integer defining how many coeff's are needed to define the diffusion
134
 
    ! COMP_DIFF_COEF( NCOMP,  NCOMP_DIFF_COEF, NPHASE  )
135
 
 
136
 
    implicit none
137
 
    INTEGER, intent( in ) :: NDIM, NPHASE, NCOMP_DIFF_COEF, &
138
 
         COMP_DIFFUSION_OPT, MAT_NONODS, TOTELE, MAT_NLOC, CV_NLOC, U_NLOC, CV_SNLOC, U_SNLOC, &
139
 
         X_NONODS, U_NONODS, CV_NONODS, U_ELE_TYPE, P_ELE_TYPE
140
 
    REAL, DIMENSION( MAT_NONODS, NDIM, NDIM, NPHASE ), intent( inout ) :: COMP_DIFFUSION
141
 
    REAL, DIMENSION( NCOMP_DIFF_COEF, NPHASE ), intent( in ) :: COMP_DIFF_COEF
142
 
    REAL, DIMENSION( X_NONODS ), intent( in ) :: X, Y, Z
143
 
    REAL, DIMENSION( U_NONODS * NPHASE ), intent( in ) :: NU, NV, NW
144
 
    INTEGER, DIMENSION( TOTELE * MAT_NLOC ), intent( in ) :: MAT_NDGLN
145
 
    INTEGER, DIMENSION( TOTELE * U_NLOC ), intent( in ) :: U_NDGLN
146
 
    INTEGER, DIMENSION( TOTELE * CV_NLOC ), intent( in ) :: X_NDGLN
147
 
 
148
 
    ! Local variables
149
 
    REAL, DIMENSION( : ), allocatable :: UD, MAT_U
150
 
    INTEGER :: ELE, CV_NOD, MAT_NOD, IPHASE, IDIM
151
 
    REAL :: DIFF_molecular, DIFF_longitudinal, DIFF_transverse
152
 
 
153
 
    IF( COMP_DIFFUSION_OPT == 0 ) THEN
154
 
       COMP_DIFFUSION = 0.0
155
 
       RETURN
156
 
    ENDIF
157
 
 
158
 
    ALLOCATE( MAT_U( NPHASE * NDIM * CV_NONODS ))
159
 
 
160
 
    CALL PROJ_U2MAT( NDIM, NPHASE, COMP_DIFFUSION_OPT, MAT_NONODS, &
161
 
         TOTELE, CV_NONODS, MAT_NLOC, CV_NLOC, U_NLOC, CV_SNLOC, U_SNLOC, &
162
 
         COMP_DIFFUSION, NCOMP_DIFF_COEF, COMP_DIFF_COEF, &
163
 
         X_NONODS, X, Y, Z, NU, NV, NW, U_NONODS, MAT_NDGLN, U_NDGLN, X_NDGLN, &
164
 
         U_ELE_TYPE, P_ELE_TYPE, &
165
 
         MAT_U ) 
166
 
 
167
 
    ! Determine the diffusion coeff tensor COMP_DIFFUSION from MAT_U and COMP_DIFF_COEF
168
 
    ALLOCATE( UD( NDIM ))
169
 
 
170
 
    DO MAT_NOD = 1, MAT_NONODS
171
 
       DO IPHASE = 1, NPHASE
172
 
 
173
 
          DO IDIM = 1,NDIM
174
 
             UD( IDIM ) = MAT_U( MAT_NOD + ( IDIM - 1 ) * MAT_NONODS + &
175
 
                  ( IPHASE - 1 ) * NDIM * MAT_NONODS )
176
 
          END DO
177
 
 
178
 
          DIFF_molecular    = COMP_DIFF_COEF( 1, IPHASE )
179
 
          DIFF_longitudinal = COMP_DIFF_COEF( 2, IPHASE )
180
 
          DIFF_transverse   = COMP_DIFF_COEF( 3, IPHASE )
181
 
 
182
 
          CALL CALC_COMP_DIF_TEN( NDIM, UD, DIFF_molecular, DIFF_longitudinal, DIFF_transverse, &
183
 
               COMP_DIFFUSION( MAT_NOD, : , : , IPHASE )) 
184
 
       END DO
185
 
    END DO
186
 
 
187
 
    DEALLOCATE( MAT_U )
188
 
    DEALLOCATE( UD )
189
 
 
190
 
    RETURN
191
 
  end subroutine CALC_COMP_DIF
192
 
 
193
 
 
194
 
 
195
 
  SUBROUTINE PROJ_U2MAT( NDIM, NPHASE, COMP_DIFFUSION_OPT, MAT_NONODS, &
196
 
       TOTELE, CV_NONODS, MAT_NLOC, CV_NLOC, U_NLOC, CV_SNLOC, U_SNLOC, &
197
 
       COMP_DIFFUSION, NCOMP_DIFF_COEF, COMP_DIFF_COEF, &
198
 
       X_NONODS, X, Y, Z, NU, NV, NW, U_NONODS, MAT_NDGLN, U_NDGLN, X_NDGLN, &
199
 
       U_ELE_TYPE, P_ELE_TYPE, &
200
 
       MAT_U ) 
201
 
    ! Determine MAT_U from NU,NV,NW which are these variables mapped to material mesh. 
202
 
    use shape_functions
203
 
    use matrix_operations
204
 
 
205
 
    implicit none
206
 
    INTEGER, intent( in ) :: NDIM, NPHASE, NCOMP_DIFF_COEF, &
207
 
         COMP_DIFFUSION_OPT, MAT_NONODS, TOTELE, MAT_NLOC, CV_NLOC, U_NLOC, CV_SNLOC, U_SNLOC, &
208
 
         X_NONODS, U_NONODS, CV_NONODS, U_ELE_TYPE, P_ELE_TYPE
209
 
    REAL, DIMENSION( MAT_NONODS, NDIM, NDIM, NPHASE ), intent( in ) :: COMP_DIFFUSION
210
 
    REAL, DIMENSION( NCOMP_DIFF_COEF, NPHASE ), intent( in ) :: COMP_DIFF_COEF
211
 
    REAL, DIMENSION( X_NONODS ), intent( in ) :: X, Y, Z
212
 
    REAL, DIMENSION( U_NONODS * NPHASE ), intent( in ) :: NU, NV, NW
213
 
    INTEGER, DIMENSION( TOTELE * MAT_NLOC ), intent( in ) :: MAT_NDGLN
214
 
    INTEGER, DIMENSION( TOTELE * U_NLOC ), intent( in ) :: U_NDGLN
215
 
    INTEGER, DIMENSION( TOTELE * CV_NLOC ), intent( in ) :: X_NDGLN
216
 
    REAL, DIMENSION( MAT_NONODS * NPHASE * NDIM), intent( inout ) :: MAT_U
217
 
    ! Determine MAT_U from NU,NV,NW which are these variables mapped to material mesh. 
218
 
 
219
 
    ! Local variables
220
 
    INTEGER, DIMENSION( :, : ), allocatable :: CV_SLOCLIST, U_SLOCLIST, CV_NEILOC, FACE_ELE
221
 
    INTEGER, DIMENSION( : ), allocatable :: FINDGPTS, COLGPTS
222
 
    REAL, DIMENSION( : ),    ALLOCATABLE :: CVWEIGHT, CVWEIGHT_SHORT, DETWEI,RA,  &
223
 
         SCVFEWEIGH, SBCVFEWEIGH, SELE_OVERLAP_SCALE
224
 
    REAL, DIMENSION( :, : ), ALLOCATABLE :: CVN, CVN_SHORT, CVFEN, CVFENLX, CVFENLY, CVFENLZ, & 
225
 
         CVFENX, CVFENY, CVFENZ, CVFEN_SHORT, CVFENLX_SHORT, CVFENLY_SHORT, CVFENLZ_SHORT, & 
226
 
         CVFENX_SHORT, CVFENY_SHORT, CVFENZ_SHORT, &
227
 
         UFEN, UFENLX, UFENLY, UFENLZ, UFENX, UFENY, UFENZ, SCVFEN, SCVFENSLX, SCVFENSLY, &
228
 
         SCVFENLX, SCVFENLY, SCVFENLZ, &
229
 
         SUFEN, SUFENSLX, SUFENSLY, SUFENLX, SUFENLY, SUFENLZ, &
230
 
         SBCVFEN, SBCVFENSLX, SBCVFENSLY, SBCVFENLX, SBCVFENLY, SBCVFENLZ, &
231
 
         SBUFEN, SBUFENSLX, SBUFENSLY, SBUFENLX, SBUFENLY, SBUFENLZ, &
232
 
         MASS, INV_MASS, MASS2U, INV_MASS_NM!, FEMU
233
 
    LOGICAL, DIMENSION( :, : ), allocatable :: CV_ON_FACE,U_ON_FACE
234
 
    INTEGER :: CV_NGI, CV_NGI_SHORT, SCVNGI, SBCVNGI, NFACE, &
235
 
         ELE, MAT_ILOC, MAT_JLOC, CV_GI, U_JLOC, MAT_KLOC, MAT_NODI, MAT_NOD, &
236
 
         U_NODJ, IPHASE, U_NODJ_IP, IDIM, JDIM, MAT_NOD_ID_IP, CV_GI_SHORT, NCOLGPTS 
237
 
    REAL :: NN, NFEMU, VOLUME, MASELE
238
 
    LOGICAL :: D1, D3, DCYL
239
 
 
240
 
 
241
 
    CALL RETRIEVE_NGI( CV_NGI, CV_NGI_SHORT, SCVNGI, SBCVNGI, NFACE, &
242
 
         NDIM, U_ELE_TYPE, CV_NLOC, U_NLOC ) 
243
 
 
244
 
    ALLOCATE( DETWEI( CV_NGI ))
245
 
    ALLOCATE( RA( CV_NGI ))
246
 
 
247
 
    ALLOCATE( CVWEIGHT( CV_NGI ))
248
 
    ALLOCATE( CVN( CV_NLOC, CV_NGI ))
249
 
    ALLOCATE( CVFEN( CV_NLOC, CV_NGI ))
250
 
    ALLOCATE( CVFENLX( CV_NLOC, CV_NGI ))
251
 
    ALLOCATE( CVFENLY( CV_NLOC, CV_NGI ))
252
 
    ALLOCATE( CVFENLZ( CV_NLOC, CV_NGI ))
253
 
    ALLOCATE( CVFENX( CV_NLOC, CV_NGI )) 
254
 
    ALLOCATE( CVFENY( CV_NLOC, CV_NGI ))
255
 
    ALLOCATE( CVFENZ( CV_NLOC, CV_NGI ))
256
 
 
257
 
    ALLOCATE( CVWEIGHT_SHORT( CV_NGI_SHORT ))
258
 
    ALLOCATE( CVN_SHORT( CV_NLOC, CV_NGI_SHORT ))
259
 
    ALLOCATE( CVFEN_SHORT( CV_NLOC, CV_NGI_SHORT))
260
 
    ALLOCATE( CVFENLX_SHORT( CV_NLOC, CV_NGI_SHORT ))
261
 
    ALLOCATE( CVFENLY_SHORT( CV_NLOC, CV_NGI_SHORT ))
262
 
    ALLOCATE( CVFENLZ_SHORT( CV_NLOC, CV_NGI_SHORT ))
263
 
    ALLOCATE( CVFENX_SHORT( CV_NLOC, CV_NGI_SHORT )) 
264
 
    ALLOCATE( CVFENY_SHORT( CV_NLOC, CV_NGI_SHORT ))
265
 
    ALLOCATE( CVFENZ_SHORT( CV_NLOC, CV_NGI_SHORT ))
266
 
 
267
 
    ALLOCATE( UFEN( U_NLOC, CV_NGI ))
268
 
    ALLOCATE( UFENLX( U_NLOC, CV_NGI ))
269
 
    ALLOCATE( UFENLY( U_NLOC, CV_NGI ))
270
 
    ALLOCATE( UFENLZ( U_NLOC, CV_NGI ))
271
 
    ALLOCATE( UFENX( U_NLOC, CV_NGI ))
272
 
    ALLOCATE( UFENY( U_NLOC, CV_NGI ))
273
 
    ALLOCATE( UFENZ( U_NLOC, CV_NGI ))
274
 
 
275
 
    ALLOCATE( SCVFEN( CV_NLOC, SCVNGI ))
276
 
    ALLOCATE( SCVFENSLX( CV_NLOC, SCVNGI ))
277
 
    ALLOCATE( SCVFENSLY( CV_NLOC, SCVNGI ))
278
 
    ALLOCATE( SCVFENLX( CV_NLOC, SCVNGI ))
279
 
    ALLOCATE( SCVFENLY( CV_NLOC, SCVNGI ))
280
 
    ALLOCATE( SCVFENLZ( CV_NLOC, SCVNGI ))
281
 
    ALLOCATE( SCVFEWEIGH( SCVNGI ))
282
 
 
283
 
    ALLOCATE( SUFEN( U_NLOC, SCVNGI ))
284
 
    ALLOCATE( SUFENSLX( U_NLOC, SCVNGI ))
285
 
    ALLOCATE( SUFENSLY( U_NLOC, SCVNGI ))
286
 
    ALLOCATE( SUFENLX( U_NLOC, SCVNGI ))
287
 
    ALLOCATE( SUFENLY( U_NLOC, SCVNGI ))
288
 
    ALLOCATE( SUFENLZ( U_NLOC, SCVNGI ))
289
 
 
290
 
    ALLOCATE( SBCVFEN( CV_SNLOC, SBCVNGI )) 
291
 
    ALLOCATE( SBCVFENSLX( CV_SNLOC, SBCVNGI ))
292
 
    ALLOCATE( SBCVFENSLY( CV_SNLOC, SBCVNGI ))
293
 
    ALLOCATE( SBCVFEWEIGH( SBCVNGI ))
294
 
    ALLOCATE( SBCVFENLX( CV_SNLOC, SBCVNGI ))
295
 
    ALLOCATE( SBCVFENLY( CV_SNLOC, SBCVNGI ))
296
 
    ALLOCATE( SBCVFENLZ( CV_SNLOC, SBCVNGI ))
297
 
    ALLOCATE( SBUFEN( U_SNLOC, SBCVNGI ))
298
 
    ALLOCATE( SBUFENSLX( U_SNLOC, SBCVNGI ))
299
 
    ALLOCATE( SBUFENSLY( U_SNLOC, SBCVNGI ))
300
 
    ALLOCATE( SBUFENLX( U_SNLOC, SBCVNGI ))
301
 
    ALLOCATE( SBUFENLY( U_SNLOC, SBCVNGI ))
302
 
    ALLOCATE( SBUFENLZ( U_SNLOC, SBCVNGI ))
303
 
 
304
 
    ALLOCATE( CV_SLOCLIST( NFACE, CV_SNLOC ))
305
 
    ALLOCATE( U_SLOCLIST( NFACE, U_SNLOC ))
306
 
    ALLOCATE( CV_NEILOC( CV_NLOC, SCVNGI ))
307
 
 
308
 
    ALLOCATE( COLGPTS( CV_NLOC * SCVNGI )) !The size of this vector is over-estimated
309
 
    ALLOCATE( FINDGPTS( CV_NLOC + 1 ))
310
 
    NCOLGPTS = 0
311
 
 
312
 
    ALLOCATE( CV_ON_FACE( CV_NLOC, SCVNGI ))
313
 
    ALLOCATE( U_ON_FACE( U_NLOC, SCVNGI ))
314
 
 
315
 
    ALLOCATE( SELE_OVERLAP_SCALE( CV_NLOC ))
316
 
 
317
 
!    ALLOCATE( FEMU( U_NLOC, CV_NGI ))
318
 
 
319
 
    CALL CV_FEM_SHAPE_FUNS( &
320
 
                                ! Volume shape functions...
321
 
         NDIM, P_ELE_TYPE,  & 
322
 
         CV_NGI, CV_NGI_SHORT, CV_NLOC, U_NLOC, CVN, CVN_SHORT, &
323
 
         CVWEIGHT, CVFEN, CVFENLX, CVFENLY, CVFENLZ, &
324
 
         CVWEIGHT_SHORT, CVFEN_SHORT, CVFENLX_SHORT, CVFENLY_SHORT, CVFENLZ_SHORT, &
325
 
         UFEN, UFENLX, UFENLY, UFENLZ, &
326
 
                                ! Surface of each CV shape functions...
327
 
         SCVNGI, CV_NEILOC, CV_ON_FACE,  &  
328
 
         SCVFEN, SCVFENSLX, SCVFENSLY, SCVFEWEIGH, &
329
 
         SCVFENLX, SCVFENLY, SCVFENLZ,  &
330
 
         SUFEN, SUFENSLX, SUFENSLY,  &
331
 
         SUFENLX, SUFENLY, SUFENLZ,  &
332
 
                                ! Surface element shape funcs...
333
 
         U_ON_FACE, NFACE, & 
334
 
         SBCVNGI, SBCVFEN, SBCVFENSLX, SBCVFENSLY, SBCVFEWEIGH, SBCVFENLX, SBCVFENLY, SBCVFENLZ, &
335
 
         SBUFEN, SBUFENSLX, SBUFENSLY, SBUFENLX, SBUFENLY, SBUFENLZ, &
336
 
         CV_SLOCLIST, U_SLOCLIST, CV_SNLOC, U_SNLOC, &
337
 
                                ! Define the gauss points that lie on the surface of the CV...
338
 
         FINDGPTS, COLGPTS, NCOLGPTS, &
339
 
         SELE_OVERLAP_SCALE ) 
340
 
 
341
 
 
342
 
    ALLOCATE( MASS( MAT_NLOC, MAT_NLOC ))
343
 
    ALLOCATE( INV_MASS( MAT_NLOC, MAT_NLOC ))
344
 
    ALLOCATE( MASS2U( MAT_NLOC, U_NLOC ))
345
 
    ALLOCATE( INV_MASS_NM( MAT_NLOC, U_NLOC ))
346
 
 
347
 
    D1 = ( NDIM == 1 )
348
 
    D3 = ( NDIM == 3 )
349
 
    DCYL = .FALSE. 
350
 
 
351
 
    Loop_Elements1: DO ELE = 1, TOTELE
352
 
 
353
 
       ! Calculate DETWEI,RA,NX,NY,NZ for element ELE
354
 
       CALL DETNLXR_PLUS_U( ELE, X, Y, Z, X_NDGLN, TOTELE, X_NONODS, CV_NLOC, CV_NGI, &
355
 
            CVFEN, CVFENLX, CVFENLY, CVFENLZ, CVWEIGHT, DETWEI, RA, VOLUME, D1, D3, DCYL, &
356
 
            CVFENX, CVFENY, CVFENZ, &
357
 
            U_NLOC, UFENLX, UFENLY, UFENLZ, UFENX, UFENY, UFENZ ) 
358
 
 
359
 
       MASELE = 0.0
360
 
       Loop_MAT_ILOC: DO MAT_ILOC = 1, MAT_NLOC
361
 
 
362
 
          Loop_MAT_JLOC: DO MAT_JLOC = 1, MAT_NLOC
363
 
 
364
 
             NN = 0.0
365
 
             DO CV_GI_SHORT = 1, CV_NGI_SHORT
366
 
                NN = NN +  CVFEN_SHORT( MAT_ILOC, CV_GI_SHORT )  * CVFEN_SHORT(  MAT_JLOC, CV_GI_SHORT ) &
367
 
                     * DETWEI( CV_GI_SHORT )
368
 
             END DO
369
 
 
370
 
             MASS( MAT_ILOC,MAT_JLOC)  = MASS( MAT_ILOC,MAT_JLOC) + NN
371
 
 
372
 
          END DO Loop_MAT_JLOC
373
 
 
374
 
       END DO Loop_MAT_ILOC
375
 
 
376
 
       CALL MATDMATINV( MASS, INV_MASS, MAT_NLOC)
377
 
 
378
 
       MASS2U = 0.0
379
 
       Loop_MAT_ILOC2: DO MAT_ILOC = 1, MAT_NLOC
380
 
 
381
 
          Loop_U_JLOC: DO U_JLOC = 1, U_NLOC
382
 
 
383
 
             NFEMU = 0.0 !!!!!!!!!!!!HHHHHHHHEEEEERRRREEEEEEEE!!!
384
 
             DO CV_GI = 1, CV_NGI
385
 
               ! NFEMU = NFEMU +  CVFEN( MAT_ILOC, CV_GI ) * FEMU(  U_JLOC, CV_GI ) * DETWEI( CV_GI )
386
 
                NFEMU = NFEMU +  CVFEN( MAT_ILOC, CV_GI ) * UFEN(  U_JLOC, CV_GI ) * DETWEI( CV_GI )
387
 
             END DO
388
 
 
389
 
             MASS2U( MAT_ILOC,U_JLOC)  = MASS2U( MAT_ILOC,U_JLOC) + NFEMU
390
 
 
391
 
          END DO Loop_U_JLOC
392
 
 
393
 
       END DO Loop_MAT_ILOC2
394
 
 
395
 
       INV_MASS_NM = 0.0
396
 
 
397
 
       DO MAT_ILOC = 1, MAT_NLOC
398
 
          DO U_JLOC = 1, U_NLOC
399
 
             DO MAT_KLOC = 1, MAT_NLOC
400
 
                INV_MASS_NM( MAT_ILOC, U_JLOC ) = INV_MASS_NM( MAT_ILOC, U_JLOC ) &
401
 
                     + INV_MASS( MAT_ILOC, MAT_KLOC ) * MASS2U( MAT_KLOC, U_JLOC )
402
 
 
403
 
             END DO
404
 
          END DO
405
 
       END DO
406
 
 
407
 
       Loop_MAT_ILOC3: DO MAT_ILOC = 1, MAT_NLOC
408
 
 
409
 
          MAT_NOD = MAT_NDGLN(( ELE - 1 ) * MAT_NLOC + MAT_ILOC )
410
 
 
411
 
          Loop_U_JLOC2: DO U_JLOC = 1, U_NLOC
412
 
 
413
 
             U_NODJ = U_NDGLN(( ELE - 1 ) * U_NLOC + U_JLOC )
414
 
 
415
 
             Loop_IPHASE: DO IPHASE = 1, NPHASE
416
 
 
417
 
                U_NODJ_IP = U_NODJ + ( IPHASE - 1 ) * U_NONODS
418
 
 
419
 
                IDIM = 1
420
 
 
421
 
                MAT_NOD_ID_IP = MAT_NOD + ( IDIM - 1 ) * MAT_NONODS  + &
422
 
                     ( IPHASE - 1 ) * NDIM * MAT_NONODS
423
 
 
424
 
                MAT_U( MAT_NOD_ID_IP ) = MAT_U( MAT_NOD_ID_IP ) + &
425
 
                     INV_MASS_NM( MAT_ILOC, U_JLOC ) * NU( U_NODJ_IP ) 
426
 
 
427
 
                IF( NDIM >= 2 ) THEN
428
 
                   IDIM = 2
429
 
                   MAT_NOD_ID_IP = MAT_NOD + ( IDIM - 1 ) * MAT_NONODS + &
430
 
                        ( IPHASE- 1 ) * NDIM * MAT_NONODS
431
 
                   MAT_U( MAT_NOD_ID_IP ) = MAT_U( MAT_NOD_ID_IP ) + &
432
 
                        INV_MASS_NM( MAT_ILOC, U_JLOC ) * NV( U_NODJ_IP ) 
433
 
                ENDIF
434
 
 
435
 
                IF( NDIM >= 3 ) THEN
436
 
                   IDIM = 3
437
 
                   MAT_NOD_ID_IP = MAT_NOD + ( IDIM - 1 ) * MAT_NONODS + &
438
 
                        ( IPHASE - 1 ) * NDIM * MAT_NONODS
439
 
                   MAT_U( MAT_NOD_ID_IP ) = MAT_U( MAT_NOD_ID_IP ) + &
440
 
                        INV_MASS_NM( MAT_ILOC, U_JLOC ) * NW( U_NODJ_IP ) 
441
 
                ENDIF
442
 
 
443
 
             END DO Loop_IPHASE
444
 
 
445
 
          END DO Loop_U_JLOC2
446
 
 
447
 
       END DO Loop_MAT_ILOC3
448
 
 
449
 
    END DO Loop_Elements1
450
 
 
451
 
 
452
 
    ! Deallocating temporary arrays
453
 
 
454
 
    DEALLOCATE( DETWEI )
455
 
    DEALLOCATE( RA )
456
 
    DEALLOCATE( CVWEIGHT )
457
 
    DEALLOCATE( CVN )
458
 
    DEALLOCATE( CVFEN )
459
 
    DEALLOCATE( CVFENLX )
460
 
    DEALLOCATE( CVFENLY )
461
 
    DEALLOCATE( CVFENLZ )
462
 
    DEALLOCATE( CVFENX ) 
463
 
    DEALLOCATE( CVFENY )
464
 
    DEALLOCATE( CVFENZ )
465
 
 
466
 
    DEALLOCATE( CVWEIGHT_SHORT )
467
 
    DEALLOCATE( CVN_SHORT )
468
 
    DEALLOCATE( CVFEN_SHORT )
469
 
    DEALLOCATE( CVFENLX_SHORT )
470
 
    DEALLOCATE( CVFENLY_SHORT )
471
 
    DEALLOCATE( CVFENLZ_SHORT )
472
 
    DEALLOCATE( CVFENX_SHORT ) 
473
 
    DEALLOCATE( CVFENY_SHORT )
474
 
    DEALLOCATE( CVFENZ_SHORT )
475
 
 
476
 
    DEALLOCATE( UFEN )
477
 
    DEALLOCATE( UFENLX )
478
 
    DEALLOCATE( UFENLY )
479
 
    DEALLOCATE( UFENLZ )
480
 
    DEALLOCATE( UFENX )
481
 
    DEALLOCATE( UFENY )
482
 
    DEALLOCATE( UFENZ )
483
 
 
484
 
    DEALLOCATE( SCVFEN )
485
 
    DEALLOCATE( SCVFENSLX )
486
 
    DEALLOCATE( SCVFENSLY )
487
 
    DEALLOCATE( SCVFENLX )
488
 
    DEALLOCATE( SCVFENLY )
489
 
    DEALLOCATE( SCVFENLZ )
490
 
    DEALLOCATE( SCVFEWEIGH )
491
 
 
492
 
    DEALLOCATE( SUFEN )
493
 
    DEALLOCATE( SUFENSLX )
494
 
    DEALLOCATE( SUFENSLY )
495
 
    DEALLOCATE( SUFENLX )
496
 
    DEALLOCATE( SUFENLY )
497
 
    DEALLOCATE( SUFENLZ )
498
 
 
499
 
    DEALLOCATE( SBCVFEN ) 
500
 
    DEALLOCATE( SBCVFENSLX )
501
 
    DEALLOCATE( SBCVFENSLY )
502
 
    DEALLOCATE( SBCVFEWEIGH )
503
 
    DEALLOCATE( SBCVFENLX )
504
 
    DEALLOCATE( SBCVFENLY )
505
 
    DEALLOCATE( SBCVFENLZ )
506
 
    DEALLOCATE( SBUFEN )
507
 
    DEALLOCATE( SBUFENSLX )
508
 
    DEALLOCATE( SBUFENSLY )
509
 
    DEALLOCATE( SBUFENLX )
510
 
    DEALLOCATE( SBUFENLY )
511
 
    DEALLOCATE( SBUFENLZ )
512
 
 
513
 
    DEALLOCATE( CV_SLOCLIST )
514
 
    DEALLOCATE( U_SLOCLIST )
515
 
    DEALLOCATE( CV_NEILOC )
516
 
 
517
 
    DEALLOCATE( COLGPTS ) !The size of this vector is over-estimated
518
 
    DEALLOCATE( FINDGPTS )
519
 
 
520
 
    DEALLOCATE( CV_ON_FACE )
521
 
    DEALLOCATE( U_ON_FACE )
522
 
 
523
 
    DEALLOCATE( SELE_OVERLAP_SCALE )
524
 
 
525
 
    DEALLOCATE( MASS )
526
 
    DEALLOCATE( INV_MASS )
527
 
    DEALLOCATE( MASS2U )
528
 
    DEALLOCATE( INV_MASS_NM )
529
 
 
530
 
    RETURN
531
 
  end subroutine PROJ_U2MAT
532
 
 
533
 
 
534
 
 
535
 
  SUBROUTINE CALC_COMP_DIF_TEN( NDIM, UD, DIFF_molecular, DIFF_longitudinal, DIFF_transverse, &
536
 
       DIFF_TEN )         
537
 
    ! Calculate the diffusion coefficient COMP_DIFFUSION for current composition...
538
 
    ! based on page 136 in Reservoir-Simulation-Mathematical-Techniques-In-Oil-Recovery-(2007).pdf
539
 
    implicit none
540
 
 
541
 
    INTEGER, intent( in ) :: NDIM
542
 
    REAL, intent( in ) :: DIFF_molecular, DIFF_longitudinal, DIFF_transverse
543
 
    REAL, DIMENSION( NDIM ), intent( in ) :: UD
544
 
    REAL, DIMENSION( NDIM, NDIM ), intent( inout ) :: DIFF_TEN
545
 
 
546
 
    ! Local variables...
547
 
    REAL, PARAMETER :: TOLER = 1.0E-10
548
 
    REAL, DIMENSION( :, : ), allocatable :: E, E_ident, E_OTH
549
 
    REAL :: RN, RN2
550
 
    INTEGER :: IDIM, JDIM
551
 
 
552
 
    ALLOCATE( E( NDIM, NDIM ))
553
 
    ALLOCATE( E_ident( NDIM, NDIM ))
554
 
    ALLOCATE( E_OTH( NDIM, NDIM ))
555
 
 
556
 
    RN = 0.0
557
 
    DO IDIM = 1, NDIM
558
 
       RN = RN + UD( IDIM ) **2
559
 
    END DO
560
 
 
561
 
    RN2 = SQRT( RN )
562
 
    RN = MAX( RN2, TOLER )
563
 
 
564
 
    E_ident = 0.0
565
 
    DO IDIM = 1, NDIM 
566
 
 
567
 
       DO JDIM = 1, NDIM
568
 
          E( IDIM, JDIM ) = UD(IDIM) * UD( JDIM ) / RN
569
 
       END DO
570
 
 
571
 
       E_ident( IDIM, IDIM ) = 1.0
572
 
    END DO
573
 
 
574
 
    E_OTH = E_ident - E
575
 
 
576
 
    DIFF_TEN = DIFF_molecular * E_ident  &
577
 
         + RN2 * ( DIFF_longitudinal * E + DIFF_transverse * E_OTH )
578
 
 
579
 
 
580
 
    DEALLOCATE( E )
581
 
    DEALLOCATE( E_ident )
582
 
    DEALLOCATE( E_OTH )
583
 
 
584
 
    RETURN
585
 
 
586
 
  END SUBROUTINE CALC_COMP_DIF_TEN
587
 
 
588
 
 
589
 
end module Compositional_Terms