APBS 3.0.0
Loading...
Searching...
No Matches
newdrvd.c
1
55#include "newdrvd.h"
56
57VEXTERNC void Vnewdriv(
58 int *iparm, double *rparm,
59 int *iwork, double *rwork,
60 double *u,
61 double *xf, double *yf, double *zf,
62 double *gxcf, double *gycf, double *gzcf,
63 double *a1cf, double *a2cf, double *a3cf,
64 double *ccf, double *fcf, double *tcf) {
65
66 int nxc;
67 int nyc;
68 int nzc;
69 int nf;
70 int nc;
71 int narr;
72 int narrc;
73 int n_rpc;
74 int n_iz;
75 int n_ipc;
76 int iretot;
77 int iintot;
78
79 int nrwk;
80 int niwk;
81 int nx;
82 int ny;
83 int nz;
84 int nlev;
85 int ierror;
86 int maxlev;
87 int mxlv;
88 int mgcoar;
89 int mgdisc;
90 int mgsolv;
91 int k_iz;
92 int k_w1;
93 int k_w2;
94 int k_ipc;
95 int k_rpc;
96 int k_ac;
97 int k_cc;
98 int k_fc;
99 int k_pc;
100
101 // Decode some parameters
102 nrwk = VAT(iparm, 1);
103 niwk = VAT(iparm, 2);
104 nx = VAT(iparm, 3);
105 ny = VAT(iparm, 4);
106 nz = VAT(iparm, 5);
107 nlev = VAT(iparm, 6);
108
109 // Some checks on input ***
110 VASSERT_MSG0(nlev > 0, "The nlev parameter must be positive");
111 VASSERT_MSG0(nx > 0, "The nx parameter must be positive");
112 VASSERT_MSG0(ny > 0, "The ny parameter must be positive");
113 VASSERT_MSG0(nz > 0, "The nz parameter must be positive");
114
115 mxlv = Vmaxlev(nx,ny,nz);
116
117 VASSERT_MSG1(nlev <= mxlv, "Max lev for your grid size is: %d", mxlv);
118
119 // Basic grid sizes, etc.
120 mgcoar = VAT(iparm, 18);
121 mgdisc = VAT(iparm, 19);
122 mgsolv = VAT(iparm, 21);
123
124 Vmgsz(&mgcoar, &mgdisc, &mgsolv,
125 &nx, &ny, &nz,
126 &nlev,
127 &nxc, &nyc, &nzc,
128 &nf, &nc,
129 &narr, &narrc,
130 &n_rpc, &n_iz, &n_ipc,
131 &iretot, &iintot);
132
133 // Allocate space for two additional work vectors ***
134 iretot = iretot + 2 * nf;
135
136 // Some more checks on input
137 VASSERT_MSG1( nrwk >= iretot, "Real work space must be: %d", iretot );
138 VASSERT_MSG1( niwk >= iintot, "Integer work space must be: %d", iintot );
139
140 // Split up the integer work array
141 k_iz = 1;
142 k_ipc = k_iz + n_iz;
143
144 // Split up the real work array
145 k_rpc = 1;
146 k_cc = k_rpc + n_rpc;
147 k_fc = k_cc + narr;
148 k_w1 = k_fc + narr;
149 k_w2 = k_w1 + nf;
150 k_pc = k_w2 + nf;
151 k_ac = k_pc + 27 * narrc;
152 // k_ac_after = 4*nf + 14*narrc;
153
154 // Call the Newton Driver
155 Vnewdriv2(iparm, rparm,
156 &nx, &ny, &nz,
157 u, RAT(iwork, k_iz),
158 RAT(rwork, k_w1), RAT(rwork, k_w2),
159 RAT(iwork, k_ipc), RAT(rwork, k_rpc),
160 RAT(rwork, k_pc), RAT(rwork, k_ac), RAT(rwork, k_cc), RAT(rwork, k_fc),
161 xf, yf, zf,
162 gxcf, gycf, gzcf,
163 a1cf, a2cf, a3cf,
164 ccf, fcf, tcf);
165}
166
167
168
169VPUBLIC void Vnewdriv2(int *iparm, double *rparm,
170 int *nx, int *ny, int *nz,
171 double *u, int *iz,
172 double *w1, double *w2,
173 int *ipc, double *rpc,
174 double *pc, double *ac, double *cc, double *fc,
175 double *xf, double *yf, double *zf,
176 double *gxcf, double *gycf, double *gzcf,
177 double *a1cf, double *a2cf, double *a3cf,
178 double *ccf, double *fcf, double *tcf) {
179
180 int mgkey;
181 int nlev;
182 int itmax;
183 int iok;
184 int iinfo;
185 int istop;
186 int ipkey;
187 int nu1;
188 int nu2;
189 int ilev;
190 int ido;
191 int iters;
192 int ierror;
193 int nlev_real;
194 int ibound;
195 int mgprol;
196 int mgcoar;
197 int mgsolv;
198 int mgdisc;
199 int mgsmoo;
200 int mode;
201 double epsiln;
202 double epsmac;
203 double errtol;
204 double omegal;
205 double omegan;
206 double bf;
207 double oh;
208 double tsetupf;
209 double tsetupc;
210 double tsolve;
211
212
213
214 // Utility variables
215 int numlev;
216
217 int iok_t;
218 int iters_t;
219 double rsnrm_t;
220 double rsden_t;
221 double orsnrm_t;
222
223 int i;
224
225
226
227 // Decode the iparm array
228 nlev = VAT(iparm, 6);
229 nu1 = VAT(iparm, 7);
230 nu2 = VAT(iparm, 8);
231 mgkey = VAT(iparm, 9);
232 itmax = VAT(iparm, 10);
233 istop = VAT(iparm, 11);
234 iinfo = VAT(iparm, 12);
235 ipkey = VAT(iparm, 14);
236 mgprol = VAT(iparm, 17);
237 mgcoar = VAT(iparm, 18);
238 mgdisc = VAT(iparm, 19);
239 mgsmoo = VAT(iparm, 20);
240 mgsolv = VAT(iparm, 21);
241
242 errtol = VAT(rparm, 1);
243 omegal = VAT(rparm, 9);
244 omegan = VAT(rparm, 10);
245
246 Vprtstp(0, -99, 0.0, 0.0, 0.0);
247
248 // Build the multigrid data structure in iz
249 Vbuildstr(nx, ny, nz, &nlev, iz);
250
251 // Start the timer
252 Vnm_tstart(30, "Vnewdrv2: fine problem setup");
253
254 // Build op and rhs on fine grid ***
255 ido = 0;
256 Vbuildops(nx, ny, nz,
257 &nlev, &ipkey, &iinfo, &ido, iz,
258 &mgprol, &mgcoar, &mgsolv, &mgdisc,
259 ipc, rpc,
260 pc, ac, cc, fc,
261 xf, yf, zf,
262 gxcf, gycf, gzcf,
263 a1cf, a2cf, a3cf,
264 ccf, fcf, tcf);
265
266 // Stop the timer
267 Vnm_tstop(30, "Vnewdrv2: fine problem setup");
268
269 // Start the timer
270 Vnm_tstart(30, "Vnewdrv2: coarse problem setup");
271
272 // Build op and rhs on all coarse grids
273 ido = 1;
274 Vbuildops(nx, ny, nz,
275 &nlev, &ipkey, &iinfo, &ido, iz,
276 &mgprol, &mgcoar, &mgsolv, &mgdisc,
277 ipc, rpc,
278 pc, ac, cc, fc,
279 xf, yf, zf,
280 gxcf, gycf, gzcf,
281 a1cf, a2cf, a3cf,
282 ccf, fcf, tcf);
283
284 // Stop the timer
285 Vnm_tstop(30, "Vnewdrv2: coarse problem setup");
286
287
288
289 /* ******************************************************************* *
290 * *** this overwrites the rhs array provided by pde specification *** *
291 * ****** compute an algebraically produced rhs for the given tcf *** *
292 * ******************************************************************* */
293 mode = 1;
294
295 if (istop == 4 || istop == 5) {
296 WARN_UNTESTED;
297 Vbuildalg(nx, ny, nz,
298 &mode, &nlev, iz,
299 ipc, rpc, ac, cc, ccf, tcf, fc, fcf);
300 }
301
302
303
304 // Determine machine epsilon
305 epsiln = Vnm_epsmac();
306
307 // Impose zero dirichlet boundary conditions (now in source fcn)
308 VfboundPMG00(nx, ny, nz, u);
309
310 // Start the timer
311 Vnm_tstart(30, "Vnewdrv2: solve");
312
313 // Call specified multigrid method
314 nlev_real = nlev;
315 iok = 1;
316 ilev = 1;
317 if (mgkey == 0) {
318 Vnewton(nx, ny, nz,
319 u, iz,
320 ccf, fcf, w1, w2,
321 &istop, &itmax, &iters, &ierror,
322 &nlev, &ilev, &nlev_real, &mgsolv,
323 &iok, &iinfo,
324 &epsiln, &errtol, &omegan,
325 &nu1, &nu2, &mgsmoo,
326 a1cf, a2cf, a3cf,
327 ipc, rpc,
328 pc, ac, cc, fc, tcf);
329 } else if (mgkey == 1) {
330 Vfnewton(nx, ny, nz,
331 u, iz, ccf, fcf, w1, w2,
332 &istop, &itmax, &iters, &ierror,
333 &nlev, &ilev, &nlev_real, &mgsolv,
334 &iok, &iinfo,
335 &epsiln, &errtol, &omegan,
336 &nu1, &nu2, &mgsmoo,
337 a1cf, a2cf, a3cf,
338 ipc, rpc,
339 pc, ac, cc, fc, tcf);
340 } else {
341 VABORT_MSG1("Bad mgkey given: %d", mgkey);
342 }
343
344 // Stop the timer
345 Vnm_tstop(30, "Vnewdrv2: solve");
346
347 // Restore boundary conditions
348 ibound = 1;
349 VfboundPMG(&ibound, nx, ny, nz, u, gxcf, gycf, gzcf);
350}
VPUBLIC void VfboundPMG(int *ibound, int *nx, int *ny, int *nz, double *x, double *gxc, double *gyc, double *gzc)
Initialize a grid function to have a certain boundary value,.
Definition mikpckd.c:209
VPUBLIC void Vbuildops(int *nx, int *ny, int *nz, int *nlev, int *ipkey, int *iinfo, int *ido, int *iz, int *mgprol, int *mgcoar, int *mgsolv, int *mgdisc, int *ipc, double *rpc, double *pc, 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, double *tcf)
Build operators, boundary arrays, modify affine vectors ido==0: do only fine level ido==1: do only co...
Definition mgsubd.c:57
VPUBLIC void Vnewton(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, int *istop, int *itmax, int *iters, int *ierror, int *nlev, int *ilev, int *nlev_real, int *mgsolv, int *iok, int *iinfo, double *epsiln, double *errtol, double *omega, int *nu1, int *nu2, int *mgsmoo, double *cprime, double *rhs, double *xtmp, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *tru)
Inexact-newton-multilevel method.
Definition newtond.c:162
VPUBLIC void Vnewdriv2(int *iparm, double *rparm, int *nx, int *ny, int *nz, double *u, int *iz, double *w1, double *w2, int *ipc, double *rpc, double *pc, 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, double *tcf)
Solves using Newton's Method.
Definition newdrvd.c:169
VEXTERNC void Vnewdriv(int *iparm, double *rparm, int *iwork, double *rwork, double *u, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf, double *tcf)
Driver for the Newton Solver.
Definition newdrvd.c:57
VPUBLIC void Vbuildstr(int *nx, int *ny, int *nz, int *nlev, int *iz)
Build the nexted operator framework in the array iz.
Definition mgsubd.c:257
VPUBLIC void Vmgsz(int *mgcoar, int *mgdisc, int *mgsolv, int *nx, int *ny, int *nz, int *nlev, int *nxc, int *nyc, int *nzc, int *nf, int *nc, int *narr, int *narrc, int *n_rpc, int *n_iz, int *n_ipc, int *iretot, int *iintot)
This routine computes the required sizes of the real and integer work arrays for the multigrid code....
Definition mgdrvd.c:562
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
Definition mikpckd.c:258
VPUBLIC void Vfnewton(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, int *istop, int *itmax, int *iters, int *ierror, int *nlev, int *ilev, int *nlev_real, int *mgsolv, int *iok, int *iinfo, double *epsiln, double *errtol, double *omega, int *nu1, int *nu2, int *mgsmoo, double *cprime, double *rhs, double *xtmp, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *tru)
Driver routines for the Newton method.
Definition newtond.c:58