APBS  3.0.0
newdrvd.c
1 
55 #include "newdrvd.h"
56 
57 VEXTERNC 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 
169 VPUBLIC 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