APBS  3.0.0
gsd.c
1 
55 #include "gsd.h"
56 
57 VPUBLIC 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 
95 VPUBLIC 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 
167 VPUBLIC 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