-
Notifications
You must be signed in to change notification settings - Fork 0
/
riccg.src
146 lines (145 loc) · 4.57 KB
/
riccg.src
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
134
135
136
137
138
139
140
141
142
143
144
145
146
#include "zeus2d.def"
c=======================================================================
c//////////////////////// SUBROUTINE RICCG \\\\\\\\\\\\\\\\\\\\\\\\\\\
c
subroutine riccg(eps,ks0,maxit)
#ifdef RAD
c
c PURPOSE:
c
c INPUT ARGUMENTS:
c eps = minimum L2 error for ICCGAF
c ks0 = max number of cyclic reductions in ICCGAF
c maxit = max number of cg iterations in ICCGAF
c
c OUTPUT ARGUMENTS:
c eps = actual L2 error achieved in ICCGAF
c maxit = actual number of cg iterations used " "
c
c EXTERNALS: ICCGAF
c
c LOCALS:
c av0,av1,bv0 = vector arrays of matrix diagonal bands. See ICCGAF
c bv1,bm1 documentation. In our application, bv1=bm1=0
c xv = vector of potential values (solution vector)
c yv = RHS of matrix eqn
c work = work space used by ICCGAF
c qa...qe = dummy variables
c-----------------------------------------------------------------------
implicit NONE
#include "param.h"
#include "grid.h"
#include "field.h"
#include "root.h"
#include "bndry.h"
#include "scratch.h"
#include "radiation.h"
REAL eps
integer maxit,ks0
integer nz
c
REAL qa,qb,qc,qd,qe
& ,av0(2*in*jn),av1(2*in*jn),bv0(2*in*jn), bv1(2*in*jn)
& ,bm1(2*in*jn),xv (2*in*jn),yv (2*in*jn),work(4*in*jn)
& ,dfn1de(in,jn),dfn1der(in,jn),dfn2de(in,jn),dfn2der(in,jn)
& ,fn1(in,jn),fn2(in,jn)
integer i,j,m
equivalence
& (av0,wa) , (av1,wc)
& , (bv0,wcg1(1 )) , (bv1 ,wcg1(2 *in*jn+1))
& , (bm1,wcg1(4*in*jn+1)) , (xv ,wcg1(6 *in*jn+1))
& , (yv ,wcg1(8*in*jn+1)) , (work,wcg1(10*in*jn+1))
c
REAL sasum
external iccgaf,sasum
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\////////////////////////////////////////
c=======================================================================
c Set up matrix elements
c
call derivs(dfn1de,dfn1der,dfn2de,dfn2der)
call rhs(fn1,fn2)
do 20 j=js,je
do 10 i=is,ie
m = (j-js)*(ie-is+1) + i - is + 1
qa = g2a(i+1)*g31a(i+1)*dx1bi(i+1)*dr1(i+1,j)
qb = g2a(i )*g31a(i )*dx1bi(i )*dr1(i ,j)
qc = g2b(i)*g2b(i)
qd = g32a(j+1)*dx2bi(j+1)*dr2(i,j+1)
qe = g32a(j )*dx2bi(j )*dr2(i,j )
c
av0(m) = -dt*(dvl1a(i)/qc*(qd+qe) + dvl2a(j)*(qa+qb))
av0(m) = av0(m) - dvl1a(i)*dvl2a(j)*(dfn1der(i,j)
& - dfn1de(i,j)*dfn2der(i,j)/dfn2de(i,j))
bv0(m) = dt*dvl1a(i)*qd/qc
av1(m) = dt*dvl2a(j)*qa
bv1(m) = 0.0
bm1(m) = 0.0
c
xv (m) = der(i,j)
yv (m) = -dvl1a(i)*dvl2a(j)
& *(dfn1de(i,j)*fn2(i,j)/dfn2de(i,j)-fn1(i,j))
10 continue
20 continue
c
c Add boundary conditions to matrix elements or yv, as appropriate.
c Order is iib, oib, ijb, ojb.
c
qa = dt*g2a(is)*g31a(is)*dx1bi(is)
do 30 j=js,je
m = (j-js)*(ie-is+1) + 1
qb = qa*dr1(is,j)
if (liib(j).eq.1 .or. liib(j).eq.2) av0(m)=av0(m)+qb*dvl2a(j)
30 continue
c
qa = dt*g2a(ie+1)*g31a(ie+1)*dx1bi(ie+1)
do 40 j=js,je
m = (j-js+1)*(ie-is+1)
av1(m) = 0.0
qb = qa*dr1(ie+1,j)
if (loib(j).eq.1 .or. loib(j).eq.2) av0(m)=av0(m)+qb*dvl2a(j)
40 continue
c
cdir$ ivdep
qa = dt*g32a(js)*dx2bi(js)
do 50 i=is,ie
m = i - is + 1
qb = dr2(i,js)*dvl1a(i)/g2b(i)**2
if (lijb(i).eq.1 .or. lijb(i).eq.2) av0(m) = av0(m) + qa*qb
50 continue
c
cdir$ ivdep
qa = dt*g32a(je+1)*dx2bi(je+1)
do 60 i=is,ie
m = (je-js)*(ie-is+1) + i - is + 1
bv0(m) = 0.0
qb = dr2(i,je+1)*dvl1a(i)/g2b(i)**2
if (lojb(i).eq.1 .or. lojb(i).eq.2) av0(m) = av0(m) + qa*qb
60 continue
c
c Matrix elements are finished, solve system.
c
c eps = eps/sasum(nx1z*nx2z,yv,1)
cRAF
cRAF Eliminate the use of pointers by aliasing some arrays through
cRAF the subroutine call -- ugly, but portable!
cRAF
cRAF call iccgaf(nx1z,nx2z,eps,ks0,maxit,av0,av1,bv0,bv1,bm1
cRAF & ,xv,yv,work)
nz = nx1z * nx2z
call iccgaf(nx1z,nx2z,eps,ks0,maxit,xv(nz+1),yv(nz+1),av0,
. work,work,work,work(nz+1),work(nz+1),
. work(2*nz+1),work(3*nz+1),
. av0,av1,bv0,bv1,bm1,xv,yv,work)
c eps = eps*sasum(nx1z*nx2z,yv,1)
c
c deconvolve xv values to get der(i,j), compute de(i,j)
c
do 80 i=is,ie
do 70 j=js,je
der(i,j) = xv((j-js)*(ie-is+1) + i - is + 1)
de (i,j) = (-fn2(i,j)-dfn2der(i,j)*der(i,j))/dfn2de(i,j)
70 continue
80 continue
#endif
return
end