← Back to branch summary

~fluidity-core/fluidity/exorcised

« back to all changes in this revision

Viewing changes to multiphase/src/multi_solvers.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_solvers.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
 
!!!=========================================!!!
29
 
!!!           SOLVERS SUBROUTINES           !!!
30
 
!!!=========================================!!!
31
 
 
32
 
 
33
 
 
34
 
module solvers_module
35
 
 
36
 
contains
37
 
 
38
 
  SUBROUTINE SOLVER( CMC, P, RHS,  &
39
 
       NCMC, NONODS, FINCMC, COLCMC, MIDCMC,  &
40
 
       ERROR, RELAX, RELAX_DIAABS, RELAX_DIA, N_LIN_ITS )
41
 
    !
42
 
    ! Solve CMC * P = RHS for RHS.
43
 
    ! RELAX: overall relaxation coeff; =1 for no relaxation. 
44
 
    ! RELAX_DIAABS: relaxation of the absolute values of the sum of the row of the matrix;
45
 
    !               - recommend >=2 for hard problems, =0 for easy
46
 
    ! RELAX_DIA: relaxation of diagonal; =1 no relaxation (normally applied). 
47
 
    ! N_LIN_ITS = no of linear iterations
48
 
    ! ERROR= solver tolerence between 2 consecutive iterations
49
 
    implicit none
50
 
    REAL, intent( in ) :: ERROR, RELAX, RELAX_DIAABS, RELAX_DIA
51
 
    INTEGER, intent( in ) ::  N_LIN_ITS, NCMC, NONODS
52
 
    REAL, DIMENSION( NCMC ), intent( in ) ::  CMC
53
 
    REAL, DIMENSION( NONODS ), intent( inout ) ::  P
54
 
    REAL, DIMENSION( NONODS ), intent( in ) :: RHS
55
 
    INTEGER, DIMENSION( NONODS + 1 ), intent( in ) :: FINCMC
56
 
    INTEGER, DIMENSION( NCMC ), intent( in ) :: COLCMC
57
 
    INTEGER, DIMENSION( NONODS ), intent( in ) :: MIDCMC
58
 
    ! Local variables
59
 
    INTEGER :: ITS, ILOOP, ISTART, IFINI, ISTEP, NOD, COUNT
60
 
    REAL :: R, SABS_DIAG, RTOP, RBOT, POLD, MAX_ERR
61
 
 
62
 
    write(357,*) 'In Solver'
63
 
 
64
 
    Loop_Non_Linear_Iter: DO ITS = 1, N_LIN_ITS
65
 
 
66
 
       MAX_ERR = 0.0
67
 
       Loop_Internal: DO ILOOP = 1, 2
68
 
          IF( ILOOP == 1 ) THEN
69
 
             ISTART = 1
70
 
             IFINI = NONODS
71
 
             ISTEP = 1
72
 
          ELSE
73
 
             ISTART = NONODS
74
 
             IFINI = 1
75
 
             ISTEP = -1
76
 
          ENDIF
77
 
 
78
 
          Loop_Nods: DO NOD = ISTART, IFINI, ISTEP
79
 
             R = RELAX_DIA * CMC( MIDCMC( NOD )) * P( NOD ) + RHS( NOD )
80
 
             SABS_DIAG = 0.0
81
 
             DO COUNT = FINCMC( NOD ), FINCMC( NOD + 1 ) - 1
82
 
                R = R - CMC( COUNT ) * P( COLCMC( COUNT ))
83
 
                SABS_DIAG = SABS_DIAG + ABS( CMC( COUNT ))
84
 
             END DO
85
 
             RTOP = R + RELAX_DIAABS * SABS_DIAG * P( NOD )
86
 
             RBOT = RELAX_DIAABS * SABS_DIAG + RELAX_DIA * CMC( MIDCMC( NOD ))
87
 
             POLD = P( NOD )
88
 
             P( NOD ) = RELAX * ( RTOP / RBOT ) + ( 1.0 - RELAX ) * P( NOD )
89
 
             MAX_ERR = MAX( MAX_ERR, ABS( POLD - P( NOD )))
90
 
          END DO Loop_Nods
91
 
       END DO Loop_Internal
92
 
 
93
 
       IF( MAX_ERR < ERROR ) CYCLE
94
 
 
95
 
    END DO Loop_Non_Linear_Iter
96
 
 
97
 
    write(357,*) 'Leaving Solver'
98
 
 
99
 
    RETURN
100
 
  END SUBROUTINE SOLVER
101
 
 
102
 
 
103
 
  SUBROUTINE GETCMC( CMC, CTP, DIAINV, CT, FREDOP, NONODS, &
104
 
       NCOLCT, FINDCT, COLCT, &
105
 
       NCMC, FINCMC, COLCMC )
106
 
    IMPLICIT NONE
107
 
    INTEGER, intent( in ) :: FREDOP, NONODS, NCOLCT, NCMC
108
 
    INTEGER, DIMENSION( NCMC ), intent( inout ) :: CMC
109
 
    INTEGER, DIMENSION( NCOLCT ), intent( in ) :: CTP
110
 
    INTEGER, DIMENSION( NONODS ), intent( in ) :: DIAINV
111
 
    INTEGER, DIMENSION( NCOLCT ), intent( in ) :: COLCT, CT
112
 
    INTEGER, DIMENSION( FREDOP + 1 ), intent( in ) :: FINDCT, FINCMC
113
 
    INTEGER, DIMENSION( NCMC ), intent( in ) :: COLCMC
114
 
 
115
 
    ! Local variables...
116
 
    INTEGER :: ROW, START, FINISH, ST2, FI2, COL, ST1, FI1, I, J, ACR, ICOL
117
 
    REAL :: CTC
118
 
 
119
 
    write(357,*) 'In GetCMC'
120
 
 
121
 
    !write(357,*)'ncmc,FREDOP,nonods=',ncmc,FREDOP,nonods
122
 
 
123
 
    CMC( 1 : NCMC ) = 0.0
124
 
 
125
 
    Loop_Row: DO ROW = 1, FREDOP
126
 
       START = FINCMC( ROW )
127
 
       FINISH = FINCMC( ROW + 1 ) - 1
128
 
       ST2 = FINDCT( ROW )
129
 
       FI2 = FINDCT( ROW + 1 ) - 1
130
 
 
131
 
       DO ACR = START, FINISH, 1
132
 
          COL = COLCMC( ACR )
133
 
          CTC = 0.
134
 
          ! Multiplying vector column 'ROW' of C1 by vector column 'COL' of C1.
135
 
          ST1 = FINDCT( COL )
136
 
          FI1 = FINDCT( COL + 1 ) - 1
137
 
 
138
 
          Loop_I: DO I = ST1, FI1
139
 
             Loop_J: DO J = ST2, FI2
140
 
                IF( COLCT( I ) == COLCT( J )) THEN
141
 
                   ICOL = COLCT( I )
142
 
                   CTC = CTC + CT( I ) * CTP( J ) * DIAINV( ICOL )
143
 
                END IF
144
 
             END DO Loop_J
145
 
          END DO Loop_I
146
 
 
147
 
          CMC( ACR ) = CTC
148
 
 
149
 
       END DO
150
 
 
151
 
    END DO Loop_Row
152
 
 
153
 
    write(357,*) 'Leaving GetCMC'
154
 
 
155
 
    RETURN
156
 
  END SUBROUTINE GETCMC
157
 
 
158
 
 
159
 
  subroutine getnc( totele, nonods, nloc, nc )
160
 
    implicit none
161
 
    integer, intent( in ) :: totele, nonods, nloc
162
 
    integer, intent( inout ) :: nc
163
 
    ! Local
164
 
    integer :: nod, ele, ele1, ele2, count
165
 
 
166
 
    write(357,*) 'In GetNC'
167
 
 
168
 
    count = 0
169
 
    Loop_cvnod: do nod = 1, nonods
170
 
       ele1 = 1 + ( nod - 2 ) / ( nloc - 1 ) 
171
 
       ele2 = 1 + ( nod - 1 ) / ( nloc - 1 ) 
172
 
 
173
 
       Loop_element: do ele = max( 1, ele1 ), min( totele, ele2 ), 1
174
 
          count = count + 1
175
 
       end do Loop_element
176
 
 
177
 
    end do Loop_cvnod
178
 
 
179
 
    nc = count
180
 
 
181
 
    write(357,*) 'Leaving GetNC'
182
 
 
183
 
    return
184
 
  end subroutine getnc
185
 
 
186
 
  subroutine getnele( band, nonods, nele )
187
 
    implicit none
188
 
    integer, intent( in ) :: band, nonods
189
 
    integer, intent( inout ) :: nele
190
 
    ! band is the semi-bandwidth. 
191
 
    ! Local
192
 
    integer :: nod, count, iband, col
193
 
 
194
 
    write(357,*) 'In GetNELE'
195
 
 
196
 
    count = 0
197
 
    col = 0
198
 
    Loop_cvnod: do nod = 1, nonods
199
 
 
200
 
       Loop_Band: do iband = -band, band, 1
201
 
          col = nod + iband
202
 
          if( ( col >= 1 ) .and. ( col <= nonods ) ) count = count + 1
203
 
       end do Loop_Band
204
 
 
205
 
    end do Loop_cvnod
206
 
 
207
 
    nele = count
208
 
 
209
 
    print *,'GetNELE'
210
 
 
211
 
    return
212
 
  end subroutine getnele
213
 
 
214
 
 
215
 
end module solvers_module
216
 
 
217