forked from shirangi/reservoir_simulator
-
Notifications
You must be signed in to change notification settings - Fork 0
/
my solver.f90
106 lines (75 loc) · 1.97 KB
/
my solver.f90
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
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!! 3D oil reservoir Simulator project
!! by Mehrdad Gharib Shirangi From Spring 2010 until Fall 2010
!! mehrdad.ghsh@gmail.com
!! Copy Right by Mehrdad Gharib Shirangi
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Subroutine solve(N,AA,BB,XX)
use globalvar
implicit none
INTEGER :: N
Real*8, intent(out) :: XX(N)
Real*8, intent(in) :: BB(N) , AA(N,N)
REAL*8, ALLOCATABLE , DIMENSION(:,:):: L , U
REAL*8, ALLOCATABLE , DIMENSION(:):: bprim
integer :: i , j , k
real*8 :: sum
! to trasform an n*n matrix into LU
Allocate (L(N,N))
Allocate (U(N,N))
Allocate (bprim(N))
L = 0
U = 0
bprim = 0
Do i = 1 ,N
L(i,1) = AA(i,1)
EndDo ! For i
!For j=1 To n
Do j = 1 ,N
U(1,j) = AA(1,j)/ L(1,1)
EndDo ! for j
!pause
! For j = 2 To n
Do j = 2 , N
!For i = j To n
Do i = j , N
sum = 0
!For k = 1 To j-1
Do k = 1 , j-1
Sum = sum + L(i,k)* U(k,j)
EndDo ! for k
L(i,j) = AA(i,j) - Sum
EndDo ! for i
U(j,j) = 1
!For i = j + 1 To n
Do i = j + 1 , N
sum = 0.0
!For k = 1 To j-1
Do k = 1 , j-1
Sum = Sum + L(j,k) * U (k,i)
EndDo ! for k
U(j,i) = (AA(j,i) - Sum) / L(j,j)
EndDo !for i
EndDo ! for j
bprim(1) = BB (1) / L(1,1)
Do i = 2 , N
Sum = 0
Do k = 1 , i-1
sum = sum + L(i,k) * bprim ( k )
End Do
bprim(i) = (BB(i) - sum) / L(i,i)
EndDo
XX(N) = bprim (N)
Do j = N - 1 , 1 , -1
sum = 0
Do k = j+1 , N
sum = sum + U(j,k) * XX(k)
EndDo
XX(j) = bprim(j) - sum
EndDo
DeAllocate (L)
DeAllocate (U)
DeAllocate (bprim)
End