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