57VPUBLIC
void Vgsrb(
int *nx,
int *ny,
int *nz,
58 int *ipc,
double *rpc,
59 double *ac,
double *cc,
double *fc,
60 double *x,
double *w1,
double *w2,
double *r,
61 int *itmax,
int *iters,
62 double *errtol,
double *omega,
63 int *iresid,
int *iadjoint) {
67 MAT2(ac, *nx * *ny * *nz, 1);
70 numdia = VAT(ipc, 11);
74 RAT2(ac, 1,1), cc, fc,
75 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
77 itmax, iters, errtol, omega, iresid, iadjoint);
78 }
else if (numdia == 27) {
81 RAT2(ac, 1, 1), cc, fc,
82 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
83 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
84 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
85 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
87 itmax, iters, errtol, omega, iresid, iadjoint);
89 Vnm_print(2,
"GSRB: invalid stencil type given...\n");
95VPUBLIC
void Vgsrb7x(
int *nx,
int *ny,
int *nz,
96 int *ipc,
double *rpc,
97 double *oC,
double *cc,
double *fc,
98 double *oE,
double *oN,
double *uC,
99 double *x,
double *w1,
double *w2,
double *r,
100 int *itmax,
int *iters,
101 double *errtol,
double *omega,
102 int *iresid,
int *iadjoint) {
106 MAT3(cc, *nx, *ny, *nz);
107 MAT3(fc, *nx, *ny, *nz);
108 MAT3( x, *nx, *ny, *nz);
109 MAT3(w1, *nx, *ny, *nz);
110 MAT3(w2, *nx, *ny, *nz);
111 MAT3( r, *nx, *ny, *nz);
113 MAT3(oE, *nx, *ny, *nz);
114 MAT3(oN, *nx, *ny, *nz);
115 MAT3(uC, *nx, *ny, *nz);
116 MAT3(oC, *nx, *ny, *nz);
118 for (*iters=1; *iters<=*itmax; (*iters)++) {
121 #pragma omp parallel for private(i, j, k, ioff)
122 for (k=2; k<=*nz-1; k++) {
123 for (j=2; j<=*ny-1; j++) {
124 ioff = (1 - *iadjoint) * ( (j + k + 2) % 2)
125 + ( *iadjoint) * (1 - (j + k + 2) % 2);
126 for (i=2+ioff; i<=*nx-1; i+=2) {
129 + VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
130 + VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
131 + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
132 + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
133 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
134 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
135 ) / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
141 #pragma omp parallel for private(i, j, k, ioff)
142 for (k=2; k<=*nz-1; k++) {
143 for (j=2; j<=*ny-1; j++) {
144 ioff = ( *iadjoint) * ( (j + k + 2) % 2 )
145 + (1 - *iadjoint) * (1 - (j + k + 2) % 2 );
146 for (i=2+ioff;i<=*nx-1; i+=2) {
149 + VAT3(oN, i, j, k) * VAT3(x, i,j+1, k)
150 + VAT3(oN, i, j-1, k) * VAT3(x, i,j-1, k)
151 + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
152 + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
153 + VAT3( uC, i, j, k-1) * VAT3(x, i, j,k-1)
154 + VAT3( uC, i, j, k) * VAT3(x, i, j,k+1)
155 ) / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
162 Vmresid7_1s(nx, ny, nz, ipc, rpc, oC, cc, fc, oE, oN, uC, x, r);
167VPUBLIC
void Vgsrb27x(
int *nx,
int *ny,
int *nz,
168 int *ipc,
double *rpc,
169 double *oC,
double *cc,
double *fc,
170 double *oE,
double *oN,
double *uC,
double *oNE,
double *oNW,
171 double *uE,
double *uW,
double *uN,
double *uS,
172 double *uNE,
double *uNW,
double *uSE,
double *uSW,
173 double *x,
double *w1,
double *w2,
double *r,
174 int *itmax,
int *iters,
175 double *errtol,
double *omega,
176 int *iresid,
int *iadjoint) {
184 double tmpO, tmpU, tmpD;
186 MAT3( cc, *nx, *ny, *nz);
187 MAT3(fc, *nx, *ny, *nz);
188 MAT3( x, *nx, *ny, *nz);
189 MAT3(w1, *nx, *ny, *nz);
190 MAT3(w2, *nx, *ny, *nz);
191 MAT3( r, *nx, *ny, *nz);
193 MAT3(oE, *nx, *ny, *nz);
194 MAT3(oN, *nx, *ny, *nz);
195 MAT3(uC, *nx, *ny, *nz);
196 MAT3(oC, *nx, *ny, *nz);
198 MAT3(oNE, *nx, *ny, *nz);
199 MAT3(oNW, *nx, *ny, *nz);
201 MAT3( uE, *nx, *ny, *nz);
202 MAT3( uW, *nx, *ny, *nz);
203 MAT3( uN, *nx, *ny, *nz);
204 MAT3( uS, *nx, *ny, *nz);
205 MAT3(uNE, *nx, *ny, *nz);
206 MAT3(uNW, *nx, *ny, *nz);
207 MAT3(uSE, *nx, *ny, *nz);
208 MAT3(uSW, *nx, *ny, *nz);
222 i1 = (1-*iadjoint) * 2 + *iadjoint * (*nx-1);
223 i2 = *iadjoint * 2 + (1-*iadjoint) * (*nx-1);
224 j1 = (1-*iadjoint) * 2 + *iadjoint * (*ny-1);
225 j2 = *iadjoint * 2 + (1-*iadjoint) * (*ny-1);
226 k1 = (1-*iadjoint) * 2 + *iadjoint * (*nz-1);
227 k2 = *iadjoint * 2 + (1-*iadjoint) * (*nz-1);
228 istep = *iadjoint*(-1) + (1-*iadjoint)*(1);
230 for (*iters=1; *iters<=*itmax; (*iters)++) {
233 for (k=2; k<=*nz-1; k++) {
235 for (j=2; j<=*ny-1; j++) {
237 ioff = (1 - *iadjoint) * ( (j + k + 2) % 2)
238 + ( *iadjoint) * (1 - (j + k + 2) % 2);
240 for (i=2+ioff; i<=*nx-1; i+=2) {
243 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
244 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
245 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
246 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
247 + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
248 + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
249 + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
250 + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
253 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
254 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
255 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
256 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
257 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
258 + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
259 + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
260 + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
261 + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
264 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
265 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
266 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
267 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
268 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
269 + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
270 + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
271 + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
272 + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
274 VAT3(x, i,j,k) = (VAT3(fc, i, j, k) + (tmpO + tmpU + tmpD))
275 / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
282 for (k=2; k<=*nz-1; k++) {
284 for (j=2; j<=*ny-1; j++) {
286 ioff = ( *iadjoint) * ( (j + k + 2) % 2)
287 + (1 - *iadjoint) * (1 - (j + k + 2) % 2);
289 for (i=2+ioff; i<=*nx-1; i+=2) {
292 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
293 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
294 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
295 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
296 + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
297 + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
298 + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
299 + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
302 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
303 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
304 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
305 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
306 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
307 + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
308 + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
309 + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
310 + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
313 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
314 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
315 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
316 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
317 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
318 + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
319 + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
320 + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
321 + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
323 VAT3(x, i,j,k) = (VAT3(fc, i, j, k) + (tmpO + tmpU + tmpD))
324 / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
332 Vmresid27_1s(nx, ny, nz,
VPUBLIC void Vgsrb(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *w1, double *w2, double *r, int *itmax, int *iters, double *errtol, double *omega, int *iresid, int *iadjoint)
Guass-Seidel solver.