1 /**
2 * @ingroup PMGC
3 * @author Tucker Beck [fortran ->c translation], Michael Holst [original]
4 * @brief
5 * @version $Id:
6 *
7 * @attention
8 * @verbatim
9 *
10 * APBS -- Adaptive Poisson-Boltzmann Solver
11 *
12 * Nathan A. Baker (nathan.baker@pnl.gov)
13 * Pacific Northwest National Laboratory
14 *
15 * Additional contributing authors listed in the code documentation.
16 *
17 * Copyright (c) 2010-2014 Battelle Memorial Institute. Developed at the Pacific Northwest National Laboratory, operated by Battelle Memorial Institute, Pacific Northwest Division for the U.S. Department Energy. Portions Copyright (c) 2002-2010, Washington University in St. Louis. Portions Copyright (c) 2002-2010, Nathan A. Baker. Portions Copyright (c) 1999-2002, The Regents of the University of California. Portions Copyright (c) 1995, Michael Holst.
18 * All rights reserved.
19 *
20 *
21 * Redistribution and use in source and binary forms, with or without
22 * modification, are permitted provided that the following conditions are met:
23 *
24 * - Redistributions of source code must retain the above copyright notice, this
25 * list of conditions and the following disclaimer.
26 *
27 * - Redistributions in binary form must reproduce the above copyright notice,
28 * this list of conditions and the following disclaimer in the documentation
29 * and/or other materials provided with the distribution.
30 *
31 * - Neither the name of Washington University in St. Louis nor the names of its
32 * contributors may be used to endorse or promote products derived from this
33 * software without specific prior written permission.
34 *
35 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
36 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
37 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
38 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
39 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
40 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
41 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
42 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
43 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
44 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
45 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
46 *
47 * @endverbatim
48 */
49
50 #include "buildAd.h"
51
VbuildA(int * nx,int * ny,int * nz,int * ipkey,int * mgdisc,int * numdia,int * ipc,double * rpc,double * ac,double * cc,double * fc,double * xf,double * yf,double * zf,double * gxcf,double * gycf,double * gzcf,double * a1cf,double * a2cf,double * a3cf,double * ccf,double * fcf)52 VPUBLIC void VbuildA(int* nx, int* ny, int* nz,
53 int* ipkey, int* mgdisc, int* numdia,
54 int* ipc, double* rpc,
55 double* ac, double* cc, double* fc,
56 double* xf, double* yf, double* zf,
57 double* gxcf, double* gycf, double* gzcf,
58 double* a1cf, double* a2cf, double* a3cf,
59 double* ccf, double* fcf) {
60
61 MAT2(ac, *nx * *ny * *nz, 14);
62
63 if (*mgdisc == 0) {
64
65 VbuildA_fv(nx, ny, nz,
66 ipkey, numdia,
67 ipc, rpc,
68 RAT2(ac, 1,1), cc, fc,
69 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
70 xf, yf, zf,
71 gxcf, gycf, gzcf,
72 a1cf, a2cf, a3cf,
73 ccf, fcf);
74
75 } else if (*mgdisc == 1) {
76
77 VbuildA_fe(nx, ny, nz,
78 ipkey, numdia,
79 ipc, rpc,
80 RAT2(ac, 1, 1), cc, fc,
81 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4), RAT2(ac, 1, 5), RAT2(ac, 1, 6),
82 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
83 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
84 xf, yf, zf,
85 gxcf, gycf, gzcf,
86 a1cf, a2cf, a3cf,
87 ccf, fcf);
88
89 } else {
90
91 Vnm_print(2, "VbuildA: Invalid discretization requested.\n");
92 /// @todo: use an APBS/MALOC method to quit/fail
93 exit(EXIT_FAILURE);
94
95 }
96 }
97
98
99
VbuildA_fv(int * nx,int * ny,int * nz,int * ipkey,int * numdia,int * ipc,double * rpc,double * oC,double * cc,double * fc,double * oE,double * oN,double * uC,double * xf,double * yf,double * zf,double * gxcf,double * gycf,double * gzcf,double * a1cf,double * a2cf,double * a3cf,double * ccf,double * fcf)100 VPUBLIC void VbuildA_fv(int *nx, int *ny, int *nz,
101 int *ipkey, int *numdia,
102 int *ipc, double *rpc,
103 double *oC, double *cc, double *fc, double *oE, double *oN, double *uC,
104 double *xf, double *yf, double *zf,
105 double *gxcf, double *gycf, double *gzcf,
106 double *a1cf, double *a2cf, double *a3cf,
107 double *ccf, double *fcf) {
108
109 int i, j, k; // @todo Document this function
110
111 /** @note: The following variables are temporaries that are not necessarily
112 * all needed as explicitly different variables. They do not
113 * persist beyond the scope of this function. This set could
114 * probably be reduced to at most 3 temporary variables. Well
115 * placed comments would go a lot farther in making the semantics
116 * of the code plain rather than producing a huge slew of seemingly
117 * homogeneous temporaries named using unclear abbreviations
118 */
119
120
121 int ike, jke, kke;
122
123
124 int nxm1, nym1, nzm1;
125
126
127 double hx, hy, hz;
128
129
130 double hxm1, hym1, hzm1;
131
132 double coef_fc;
133
134 double bc_cond_e;
135 double bc_cond_w;
136 double bc_cond_n;
137 double bc_cond_s;
138 double bc_cond_u;
139 double bc_cond_d;
140 double coef_oE;
141 double coef_oN;
142 double coef_uC;
143 double coef_oEm1;
144 double coef_oNm1;
145 double coef_uCm1;
146
147 double diag;
148
149 MAT3( fc, *nx, *ny, *nz);
150 MAT3( fcf, *nx, *ny, *nz);
151 MAT3( cc, *nx, *ny, *nz);
152 MAT3( ccf, *nx, *ny, *nz);
153 MAT3( oC, *nx, *ny, *nz);
154 MAT3(a1cf, *nx, *ny, *nz);
155 MAT3(a2cf, *nx, *ny, *nz);
156 MAT3(a3cf, *nx, *ny, *nz);
157 MAT3( uC, *nx, *ny, *nz);
158 MAT3( oN, *nx, *ny, *nz);
159 MAT3( oE, *nx, *ny, *nz);
160 MAT3(gxcf, *ny, *nz, 2);
161 MAT3(gycf, *nx, *nz, 2);
162 MAT3(gzcf, *nx, *ny, 2);
163
164 // Save the problem key with this operator. @todo: What?
165 VAT(ipc, 10) = *ipkey;
166
167 // Note how many nonzeros in this discretization stencil
168 VAT(ipc, 11) = 7;
169 VAT(ipc, 12) = 1;
170 *numdia = 4;
171
172 // Define n and determine number of mesh points
173 nxm1 = *nx - 1;
174 nym1 = *ny - 1;
175 nzm1 = *nz - 1;
176
177 // Determine diag scale factor
178 // (would like something close to ones on the main diagonal)
179 // @todo: Make a more meaningful comment
180 diag = 1.0;
181
182
183
184 /* *********************************************************************
185 * *** interior points ***
186 * ********************************************************************* */
187
188 // build the operator
189 //fprintf(data, "%s\n", PRINT_FUNC);
190 for(k=2; k<=*nz-1; k++) {
191
192 hzm1 = VAT(zf, k) - VAT(zf, k-1);
193 hz = VAT(zf, k+1) - VAT(zf, k);
194
195 for(j=2; j<=*ny-1; j++) {
196
197 hym1 = VAT(yf, j) - VAT(yf, j-1);
198 hy = VAT(yf, j+1) - VAT(yf, j);
199
200 for(i=2; i<=*nx-1; i++) {
201
202 hxm1 = VAT(xf, i) - VAT(xf, i-1);
203 hx = VAT(xf, i+1) - VAT(xf, i);
204
205 // Calculate some coefficients
206 /** @note that these and the running OC calculation could
207 * easily be pushed down into the step by step
208 * compuation of neighbors. That would alleviate some
209 * of the temporary madness
210 */
211
212 coef_oE = diag * (hym1 + hy) * (hzm1 + hz) / (4.0 * hx);
213 coef_oEm1 = diag * (hym1 + hy) * (hzm1 + hz) / (4.0 * hxm1);
214 coef_oN = diag * (hxm1 + hx) * (hzm1 + hz) / (4.0 * hy);
215 coef_oNm1 = diag * (hxm1 + hx) * (hzm1 + hz) / (4.0 * hym1);
216 coef_uC = diag * (hxm1 + hx) * (hym1 + hy) / (4.0 * hz);
217 coef_uCm1 = diag * (hxm1 + hx) * (hym1 + hy) / (4.0 * hzm1);
218 coef_fc = diag * (hxm1 + hx) * (hym1 + hy) * (hzm1 + hz) / 8.0;
219
220 // Calculate the coefficient and source function
221 VAT3(fc, i, j, k) = coef_fc * VAT3(fcf, i, j, k);
222 VAT3(cc, i, j, k) = coef_fc * VAT3(ccf, i, j, k);
223 //fprintf(data, "%19.12E\n", VAT3(cc, i, j, k));
224
225 // Calculate the diagonal for matvecs and smoothings
226
227 VAT3(oC, i, j, k) = coef_oE * VAT3(a1cf, i, j, k) +
228 coef_oEm1 * VAT3(a1cf, i-1, j, k) +
229 coef_oN * VAT3(a2cf, i, j, k) +
230 coef_oNm1 * VAT3(a2cf, i, j-1, k) +
231 coef_uC * VAT3(a3cf, i, j, k) +
232 coef_uCm1 * VAT3(a3cf, i, j, k-1);
233
234 //fprintf(data, "%19.12E\n", VAT3(oC, i, j, k));
235
236 // Calculate the east neighbor
237 ike = VMIN2(1, VABS(i - nxm1));
238 VAT3(oE, i, j, k) = ike * coef_oE * VAT3(a1cf, i, j, k);
239 //fprintf(data, "%19.12E\n", VAT3(oE, i, j, k));
240 bc_cond_e = (1 - ike) * coef_oE * VAT3(a1cf, i, j, k) * VAT3(gxcf, j, k, 2);
241 VAT3(fc, i, j, k) += bc_cond_e;
242
243 // Calculate the north neighbor
244 jke = VMIN2(1, VABS(j - nym1));
245 VAT3(oN, i, j, k) = jke * coef_oN * VAT3(a2cf, i, j, k);
246 //fprintf(data, "%19.12E\n", VAT3(oN, i, j, k));
247 bc_cond_n = (1 - jke) * coef_oN * VAT3(a2cf, i, j, k) * VAT3(gycf, i, k, 2);
248 VAT3(fc, i, j, k) += bc_cond_n;
249
250 // Calculate the up neighbor
251 kke = VMIN2(1, VABS(k - nzm1));
252 VAT3(uC, i, j, k) = kke * coef_uC * VAT3(a3cf, i, j, k);
253 //fprintf(data, "%19.12E\n", VAT3(uC, i, j, k));
254 bc_cond_u = (1 - kke) * coef_uC * VAT3(a3cf, i, j, k) * VAT3(gzcf, i, j, 2);
255 VAT3(fc, i, j, k) += bc_cond_u;
256
257 // Calculate the west neighbor (just handle b.c.)
258 ike = VMIN2(1, VABS(i - 2));
259 bc_cond_w = (1 - ike) * coef_oEm1 * VAT3(a1cf, i-1, j, k) * VAT3(gxcf, j, k, 1);
260 VAT3(fc, i, j, k) += bc_cond_w;
261
262 // Calculate the south neighbor (just handle b.c.)
263 jke = VMIN2(1, VABS(j - 2));
264 bc_cond_s = (1 - jke) * coef_oNm1 * VAT3(a2cf, i, j-1, k) * VAT3(gycf, i, k, 1);
265 VAT3(fc, i, j, k) += bc_cond_s;
266
267 // Calculate the down neighbor (just handle b.c.)
268 kke = VMIN2(1, VABS(k - 2));
269 bc_cond_d = (1 - kke) * coef_uCm1 * VAT3(a3cf, i, j, k-1) * VAT3(gzcf, i, j, 1);
270 VAT3(fc, i, j, k) += bc_cond_d;
271
272 //fprintf(data, "%19.12E\n", VAT3(fc, i, j, k));
273 }
274 }
275 }
276 }
277
278
279
VbuildA_fe(int * nx,int * ny,int * nz,int * ipkey,int * numdia,int * ipc,double * rpc,double * oC,double * cc,double * fc,double * oE,double * oN,double * uC,double * oNE,double * oNW,double * uE,double * uW,double * uN,double * uS,double * uNE,double * uNW,double * uSE,double * uSW,double * xf,double * yf,double * zf,double * gxcf,double * gycf,double * gzcf,double * a1cf,double * a2cf,double * a3cf,double * ccf,double * fcf)280 VPUBLIC void VbuildA_fe(int *nx, int *ny, int *nz,
281 int *ipkey, int *numdia,
282 int *ipc, double *rpc,
283 double *oC, double *cc, double *fc,
284 double *oE, double *oN, double *uC,
285 double* oNE, double* oNW,
286 double* uE, double* uW,
287 double* uN, double* uS,
288 double* uNE, double* uNW, double* uSE, double* uSW,
289 double *xf, double *yf, double *zf,
290 double *gxcf, double *gycf, double *gzcf,
291 double *a1cf, double *a2cf, double *a3cf,
292 double *ccf, double *fcf) {
293 VABORT_MSG0("Untranslated Component: from buildAd.f");
294 }
295