APBS 3.0.0
Loading...
Searching...
No Matches
gsd.c
1
55#include "gsd.h"
56
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) {
64
65 int numdia;
66
67 MAT2(ac, *nx * *ny * *nz, 1);
68
69 // Do in one step ***
70 numdia = VAT(ipc, 11);
71 if (numdia == 7) {
72 Vgsrb7x(nx, ny, nz,
73 ipc, rpc,
74 RAT2(ac, 1,1), cc, fc,
75 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
76 x, w1, w2, r,
77 itmax, iters, errtol, omega, iresid, iadjoint);
78 } else if (numdia == 27) {
79 Vgsrb27x(nx, ny, nz,
80 ipc, rpc,
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),
86 x, w1, w2, r,
87 itmax, iters, errtol, omega, iresid, iadjoint);
88 } else {
89 Vnm_print(2, "GSRB: invalid stencil type given...\n");
90 }
91}
92
93
94
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) {
103
104 int i, j, k, ioff;
105
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);
112
113 MAT3(oE, *nx, *ny, *nz);
114 MAT3(oN, *nx, *ny, *nz);
115 MAT3(uC, *nx, *ny, *nz);
116 MAT3(oC, *nx, *ny, *nz);
117
118 for (*iters=1; *iters<=*itmax; (*iters)++) {
119
120 // Do the red points ***
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) {
127 VAT3(x, i, j, k) = (
128 VAT3(fc, i, j, k)
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));
136 }
137 }
138 }
139
140 // Do the black points
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) {
147 VAT3(x, i, j, k) = (
148 VAT3(fc, i, j, k)
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));
156 }
157 }
158 }
159 }
160
161 if (*iresid == 1)
162 Vmresid7_1s(nx, ny, nz, ipc, rpc, oC, cc, fc, oE, oN, uC, x, r);
163}
164
165
166
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) {
177
178 int i, j, k;
179 int i1, j1, k1;
180 int i2, j2, k2;
181 int ioff;
182 int istep;
183
184 double tmpO, tmpU, tmpD;
185
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);
192
193 MAT3(oE, *nx, *ny, *nz);
194 MAT3(oN, *nx, *ny, *nz);
195 MAT3(uC, *nx, *ny, *nz);
196 MAT3(oC, *nx, *ny, *nz);
197
198 MAT3(oNE, *nx, *ny, *nz);
199 MAT3(oNW, *nx, *ny, *nz);
200
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);
209
210 // Do the gauss-seidel iteration itmax times
211
212 /*
213 i1 = (1 - *iadjoint) * 2 + ( *iadjoint) * (*nx - 1);
214 i2 = ( *iadjoint) * 2 + (1 - *iadjoint) * (*nx - 1);
215 j1 = (1 - *iadjoint) * 2 + ( *iadjoint) * (*ny - 1);
216 j2 = ( *iadjoint) * 2 + (1 - *iadjoint) * (*ny - 1);
217 k1 = (1 - *iadjoint) * 2 + ( *iadjoint) * (*nz - 1);
218 k2 = ( *iadjoint) * 2 + (1 - *iadjoint) * (*nz - 1);
219 istep = ( *iadjoint) * (-1) + (1 - *iadjoint) * (1);
220 */
221
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);
229
230 for (*iters=1; *iters<=*itmax; (*iters)++) {
231
232 //#pragma omp parallel for private(i, j, k, ioff, tmpO, tmpU, tmpD)
233 for (k=2; k<=*nz-1; k++) {
234
235 for (j=2; j<=*ny-1; j++) {
236
237 ioff = (1 - *iadjoint) * ( (j + k + 2) % 2)
238 + ( *iadjoint) * (1 - (j + k + 2) % 2);
239
240 for (i=2+ioff; i<=*nx-1; i+=2) {
241
242 tmpO =
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);
251
252 tmpU =
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);
262
263 tmpD =
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);
273
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));
276
277 }
278 }
279 }
280
281 //#pragma omp parallel for private(i, j, k, ioff, tmpO, tmpU, tmpD)
282 for (k=2; k<=*nz-1; k++) {
283
284 for (j=2; j<=*ny-1; j++) {
285
286 ioff = ( *iadjoint) * ( (j + k + 2) % 2)
287 + (1 - *iadjoint) * (1 - (j + k + 2) % 2);
288
289 for (i=2+ioff; i<=*nx-1; i+=2) {
290
291 tmpO =
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);
300
301 tmpU =
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);
311
312 tmpD =
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);
322
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));
325 }
326 }
327 }
328 }
329
330 // If specified, return the new residual as well
331 if (*iresid == 1)
332 Vmresid27_1s(nx, ny, nz,
333 ipc, rpc,
334 oC, cc, fc,
335 oE, oN, uC,
336 oNE, oNW,
337 uE, uW, uN, uS,
338 uNE, uNW, uSE, uSW,
339 x, r);
340}
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.
Definition gsd.c:57