APBS  1.5
gsd.c
1 
50 #include "gsd.h"
51 
52 VPUBLIC void Vgsrb(int *nx, int *ny, int *nz,
53  int *ipc, double *rpc,
54  double *ac, double *cc, double *fc,
55  double *x, double *w1, double *w2, double *r,
56  int *itmax, int *iters,
57  double *errtol, double *omega,
58  int *iresid, int *iadjoint) {
59 
60  int numdia;
61 
62  MAT2(ac, *nx * *ny * *nz, 1);
63 
64  // Do in one step ***
65  numdia = VAT(ipc, 11);
66  if (numdia == 7) {
67  Vgsrb7x(nx, ny, nz,
68  ipc, rpc,
69  RAT2(ac, 1,1), cc, fc,
70  RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
71  x, w1, w2, r,
72  itmax, iters, errtol, omega, iresid, iadjoint);
73  } else if (numdia == 27) {
74  Vgsrb27x(nx, ny, nz,
75  ipc, rpc,
76  RAT2(ac, 1, 1), cc, fc,
77  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
78  RAT2(ac, 1, 5), RAT2(ac, 1, 6),
79  RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
80  RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
81  x, w1, w2, r,
82  itmax, iters, errtol, omega, iresid, iadjoint);
83  } else {
84  Vnm_print(2, "GSRB: invalid stencil type given...\n");
85  }
86 }
87 
88 
89 
90 VPUBLIC void Vgsrb7x(int *nx,int *ny,int *nz,
91  int *ipc, double *rpc,
92  double *oC, double *cc, double *fc,
93  double *oE, double *oN, double *uC,
94  double *x, double *w1, double *w2, double *r,
95  int *itmax, int *iters,
96  double *errtol, double *omega,
97  int *iresid, int *iadjoint) {
98 
99  int i, j, k, ioff;
100 
101  MAT3(cc, *nx, *ny, *nz);
102  MAT3(fc, *nx, *ny, *nz);
103  MAT3( x, *nx, *ny, *nz);
104  MAT3(w1, *nx, *ny, *nz);
105  MAT3(w2, *nx, *ny, *nz);
106  MAT3( r, *nx, *ny, *nz);
107 
108  MAT3(oE, *nx, *ny, *nz);
109  MAT3(oN, *nx, *ny, *nz);
110  MAT3(uC, *nx, *ny, *nz);
111  MAT3(oC, *nx, *ny, *nz);
112 
113  for (*iters=1; *iters<=*itmax; (*iters)++) {
114 
115  // Do the red points ***
116  #pragma omp parallel for private(i, j, k, ioff)
117  for (k=2; k<=*nz-1; k++) {
118  for (j=2; j<=*ny-1; j++) {
119  ioff = (1 - *iadjoint) * ( (j + k + 2) % 2)
120  + ( *iadjoint) * (1 - (j + k + 2) % 2);
121  for (i=2+ioff; i<=*nx-1; i+=2) {
122  VAT3(x, i, j, k) = (
123  VAT3(fc, i, j, k)
124  + VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
125  + VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
126  + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
127  + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
128  + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
129  + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
130  ) / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
131  }
132  }
133  }
134 
135  // Do the black points
136  #pragma omp parallel for private(i, j, k, ioff)
137  for (k=2; k<=*nz-1; k++) {
138  for (j=2; j<=*ny-1; j++) {
139  ioff = ( *iadjoint) * ( (j + k + 2) % 2 )
140  + (1 - *iadjoint) * (1 - (j + k + 2) % 2 );
141  for (i=2+ioff;i<=*nx-1; i+=2) {
142  VAT3(x, i, j, k) = (
143  VAT3(fc, i, j, k)
144  + VAT3(oN, i, j, k) * VAT3(x, i,j+1, k)
145  + VAT3(oN, i, j-1, k) * VAT3(x, i,j-1, k)
146  + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
147  + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
148  + VAT3( uC, i, j, k-1) * VAT3(x, i, j,k-1)
149  + VAT3( uC, i, j, k) * VAT3(x, i, j,k+1)
150  ) / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
151  }
152  }
153  }
154  }
155 
156  if (*iresid == 1)
157  Vmresid7_1s(nx, ny, nz, ipc, rpc, oC, cc, fc, oE, oN, uC, x, r);
158 }
159 
160 
161 
162 VPUBLIC void Vgsrb27x(int *nx,int *ny,int *nz,
163  int *ipc, double *rpc,
164  double *oC, double *cc, double *fc,
165  double *oE, double *oN, double *uC, double *oNE, double *oNW,
166  double *uE, double *uW, double *uN, double *uS,
167  double *uNE, double *uNW, double *uSE, double *uSW,
168  double *x, double *w1, double *w2, double *r,
169  int *itmax, int *iters,
170  double *errtol, double *omega,
171  int *iresid, int *iadjoint) {
172 
173  int i, j, k;
174  int i1, j1, k1;
175  int i2, j2, k2;
176  int ioff;
177  int istep;
178 
179  double tmpO, tmpU, tmpD;
180 
181  MAT3( cc, *nx, *ny, *nz);
182  MAT3(fc, *nx, *ny, *nz);
183  MAT3( x, *nx, *ny, *nz);
184  MAT3(w1, *nx, *ny, *nz);
185  MAT3(w2, *nx, *ny, *nz);
186  MAT3( r, *nx, *ny, *nz);
187 
188  MAT3(oE, *nx, *ny, *nz);
189  MAT3(oN, *nx, *ny, *nz);
190  MAT3(uC, *nx, *ny, *nz);
191  MAT3(oC, *nx, *ny, *nz);
192 
193  MAT3(oNE, *nx, *ny, *nz);
194  MAT3(oNW, *nx, *ny, *nz);
195 
196  MAT3( uE, *nx, *ny, *nz);
197  MAT3( uW, *nx, *ny, *nz);
198  MAT3( uN, *nx, *ny, *nz);
199  MAT3( uS, *nx, *ny, *nz);
200  MAT3(uNE, *nx, *ny, *nz);
201  MAT3(uNW, *nx, *ny, *nz);
202  MAT3(uSE, *nx, *ny, *nz);
203  MAT3(uSW, *nx, *ny, *nz);
204 
205  // Do the gauss-seidel iteration itmax times
206 
207  /*
208  i1 = (1 - *iadjoint) * 2 + ( *iadjoint) * (*nx - 1);
209  i2 = ( *iadjoint) * 2 + (1 - *iadjoint) * (*nx - 1);
210  j1 = (1 - *iadjoint) * 2 + ( *iadjoint) * (*ny - 1);
211  j2 = ( *iadjoint) * 2 + (1 - *iadjoint) * (*ny - 1);
212  k1 = (1 - *iadjoint) * 2 + ( *iadjoint) * (*nz - 1);
213  k2 = ( *iadjoint) * 2 + (1 - *iadjoint) * (*nz - 1);
214  istep = ( *iadjoint) * (-1) + (1 - *iadjoint) * (1);
215  */
216 
217  i1 = (1-*iadjoint) * 2 + *iadjoint * (*nx-1);
218  i2 = *iadjoint * 2 + (1-*iadjoint) * (*nx-1);
219  j1 = (1-*iadjoint) * 2 + *iadjoint * (*ny-1);
220  j2 = *iadjoint * 2 + (1-*iadjoint) * (*ny-1);
221  k1 = (1-*iadjoint) * 2 + *iadjoint * (*nz-1);
222  k2 = *iadjoint * 2 + (1-*iadjoint) * (*nz-1);
223  istep = *iadjoint*(-1) + (1-*iadjoint)*(1);
224 
225  for (*iters=1; *iters<=*itmax; (*iters)++) {
226 
227  //#pragma omp parallel for private(i, j, k, ioff, tmpO, tmpU, tmpD)
228  for (k=2; k<=*nz-1; k++) {
229 
230  for (j=2; j<=*ny-1; j++) {
231 
232  ioff = (1 - *iadjoint) * ( (j + k + 2) % 2)
233  + ( *iadjoint) * (1 - (j + k + 2) % 2);
234 
235  for (i=2+ioff; i<=*nx-1; i+=2) {
236 
237  tmpO =
238  + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
239  + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
240  + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
241  + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
242  + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
243  + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
244  + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
245  + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
246 
247  tmpU =
248  + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
249  + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
250  + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
251  + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
252  + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
253  + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
254  + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
255  + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
256  + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
257 
258  tmpD =
259  + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
260  + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
261  + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
262  + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
263  + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
264  + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
265  + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
266  + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
267  + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
268 
269  VAT3(x, i,j,k) = (VAT3(fc, i, j, k) + (tmpO + tmpU + tmpD))
270  / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
271 
272  }
273  }
274  }
275 
276  //#pragma omp parallel for private(i, j, k, ioff, tmpO, tmpU, tmpD)
277  for (k=2; k<=*nz-1; k++) {
278 
279  for (j=2; j<=*ny-1; j++) {
280 
281  ioff = ( *iadjoint) * ( (j + k + 2) % 2)
282  + (1 - *iadjoint) * (1 - (j + k + 2) % 2);
283 
284  for (i=2+ioff; i<=*nx-1; i+=2) {
285 
286  tmpO =
287  + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
288  + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
289  + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
290  + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
291  + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
292  + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
293  + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
294  + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
295 
296  tmpU =
297  + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
298  + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
299  + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
300  + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
301  + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
302  + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
303  + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
304  + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
305  + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
306 
307  tmpD =
308  + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
309  + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
310  + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
311  + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
312  + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
313  + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
314  + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
315  + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
316  + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
317 
318  VAT3(x, i,j,k) = (VAT3(fc, i, j, k) + (tmpO + tmpU + tmpD))
319  / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
320  }
321  }
322  }
323  }
324 
325  // If specified, return the new residual as well
326  if (*iresid == 1)
327  Vmresid27_1s(nx, ny, nz,
328  ipc, rpc,
329  oC, cc, fc,
330  oE, oN, uC,
331  oNE, oNW,
332  uE, uW, uN, uS,
333  uNE, uNW, uSE, uSW,
334  x, r);
335 }
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:52
int ny
Definition: vpmgp.h:84
int nx
Definition: vpmgp.h:83
int itmax
Definition: vpmgp.h:122
double errtol
Definition: vpmgp.h:121
int nz
Definition: vpmgp.h:85