APBS  1.4.1
vpmgp.c
Go to the documentation of this file.
1 
57 #include "vpmgp.h"
58 
59 VEMBED(rcsid="$Id$")
60 
61 /* ///////////////////////////////////////////////////////////////////////////
62 // Class Vpmgp: Inlineable methods
64 #if !defined(VINLINE_VACC)
65 #endif /* if !defined(VINLINE_VACC) */
66 
67 /* ///////////////////////////////////////////////////////////////////////////
68 // Class Vpmgp: Non-inlineable methods
70 
71 /* ///////////////////////////////////////////////////////////////////////////
72 // Routine: Vpmgp_ctor
73 //
74 // Author: Nathan Baker
76 VPUBLIC Vpmgp* Vpmgp_ctor(MGparm *mgparm) {
77 
78  Vpmgp *thee = VNULL;
79 
80  /* Set up the structure */
81  thee = (Vpmgp*)Vmem_malloc(VNULL, 1, sizeof(Vpmgp) );
82  VASSERT( thee != VNULL);
83  VASSERT(Vpmgp_ctor2(thee,mgparm));
84 
85  return thee;
86 }
87 
88 /* ///////////////////////////////////////////////////////////////////////////
89 // Routine: Vpmgp_ctor2
90 //
91 // Author: Nathan Baker
93 VPUBLIC int Vpmgp_ctor2(Vpmgp *thee,MGparm *mgparm) {
94 
95  /* Specified parameters */
96  thee->nx = mgparm->dime[0];
97  thee->ny = mgparm->dime[1];
98  thee->nz = mgparm->dime[2];
99  thee->hx = mgparm->grid[0];
100  thee->hy = mgparm->grid[1];
101  thee->hzed = mgparm->grid[2];
102  thee->xlen = ((double)(mgparm->dime[0]-1))*mgparm->grid[0];
103  thee->ylen = ((double)(mgparm->dime[1]-1))*mgparm->grid[1];
104  thee->zlen = ((double)(mgparm->dime[2]-1))*mgparm->grid[2];
105  thee->nlev = mgparm->nlev;
106 
107  thee->nonlin = mgparm->nonlintype;
108  thee->meth = mgparm->method;
109 
110 #ifdef DEBUG_MAC_OSX_OCL
111 #include "mach_chud.h"
112  if(kOpenCLAvailable)
113  thee->meth = 4;
114 #endif
115 
116  if (thee->nonlin == NONLIN_LPBE) thee->ipkey = IPKEY_LPBE; /* LPBE case */
117  else if(thee->nonlin == NONLIN_SMPBE) thee->ipkey = IPKEY_SMPBE; /* SMPBE case */
118  else thee->ipkey = IPKEY_NPBE; /* NPBE standard case */
119 
120  /* Default parameters */
121  if (mgparm->setetol) { /* If etol is set by the user in APBS input file, then use this custom-defined etol */
122  thee->errtol = mgparm->etol;
123  Vnm_print(1, " Error tolerance (etol) is now set to user-defined \
124 value: %g \n", thee->errtol);
125  Vnm_print(0, "Error tolerance (etol) is now set to user-defined \
126 value: %g \n", thee->errtol);
127  } else thee->errtol = 1.0e-6; /* Here are a few comments. Mike had this set to
128  * 1e-9; convential wisdom sets this at 1e-6 for
129  * the PBE; Ray Luo sets this at 1e-3 for his
130  * accelerated PBE (for dynamics, etc.) */
131  thee->itmax = 200;
132  thee->istop = 1;
133  thee->iinfo = 1; /* I'd recommend either 1 (for debugging LPBE) or 2 (for debugging NPBE), higher values give too much output */
134 
135  thee->bcfl = BCFL_SDH;
136  thee->key = 0;
137  thee->iperf = 0;
138  thee->mgcoar = 2;
139  thee->mgkey = 0;
140  thee->nu1 = 2;
141  thee->nu2 = 2;
142  thee->mgprol = 0;
143  thee->mgdisc = 0;
144  thee->omegal = 19.4e-1;
145  thee->omegan = 9.0e-1;
146  thee->ipcon = 3;
147  thee->irite = 8;
148  thee->xcent = 0.0;
149  thee->ycent = 0.0;
150  thee->zcent = 0.0;
151 
152  /* Default value for all APBS runs */
153  thee->mgsmoo = 1;
154  if (thee->nonlin == NONLIN_NPBE || thee->nonlin == NONLIN_SMPBE) {
155  /* SMPBE Added - SMPBE needs to mimic NPBE */
156  Vnm_print(0, "Vpmp_ctor2: Using meth = 1, mgsolv = 0\n");
157  thee->mgsolv = 0;
158  } else {
159  /* Most rigorous (good for testing) */
160  Vnm_print(0, "Vpmp_ctor2: Using meth = 2, mgsolv = 1\n");
161  thee->mgsolv = 1;
162  }
163 
164  /* TEMPORARY USEAQUA */
165  /* If we are using aqua, our solution method is either VSOL_CGMGAqua or VSOL_NewtonAqua
166  * so we need to temporarily override the mgsolve value and set it to 0
167  */
168  if(mgparm->useAqua == 1) thee->mgsolv = 0;
169 
170  return 1;
171 }
172 
173 /* ///////////////////////////////////////////////////////////////////////////
174 // Routine: Vpmgp_dtor
175 //
176 // Author: Nathan Baker
178 VPUBLIC void Vpmgp_dtor(Vpmgp **thee) {
179 
180  if ((*thee) != VNULL) {
181  Vpmgp_dtor2(*thee);
182  Vmem_free(VNULL, 1, sizeof(Vpmgp), (void **)thee);
183  (*thee) = VNULL;
184  }
185 
186 }
187 
188 /* ///////////////////////////////////////////////////////////////////////////
189 // Routine: Vpmgp_dtor2
190 //
191 // Author: Nathan Baker
193 VPUBLIC void Vpmgp_dtor2(Vpmgp *thee) { ; }
194 
195 
196 VPUBLIC void Vpmgp_size(
197  Vpmgp *thee
198  )
199 {
200 
201  int num_nf = 0;
202  int num_narr = 2;
203  int num_narrc = 27;
204  int nxf, nyf, nzf, level, num_nf_oper, num_narrc_oper, n_band, nc_band, num_band, iretot;
205 
206  thee->nf = thee->nx * thee->ny * thee->nz;
207  thee->narr = thee->nf;
208  nxf = thee->nx;
209  nyf = thee->ny;
210  nzf = thee->nz;
211  thee->nxc = thee->nx;
212  thee->nyc = thee->ny;
213  thee->nzc = thee->nz;
214 
215  for (level=2; level<=thee->nlev; level++) {
216  Vpmgp_makeCoarse(1, nxf, nyf, nzf, &(thee->nxc), &(thee->nyc), &(thee->nzc)); /* NAB TO-DO -- implement this function and check which variables need to be passed by reference... */
217  nxf = thee->nxc;
218  nyf = thee->nyc;
219  nzf = thee->nzc;
220  thee->narr = thee->narr + (nxf * nyf * nzf);
221  }
222 
223  thee->nc = thee->nxc * thee->nyc * thee->nzc;
224  thee->narrc = thee->narr - thee->nf;
225 
226  /* Box or FEM discretization on fine grid? */
227  switch (thee->mgdisc) { /* NAB TO-DO: This needs to be changed into an enumeration */
228  case 0:
229  num_nf_oper = 4;
230  break;
231  case 1:
232  num_nf_oper = 14;
233  break;
234  default:
235  Vnm_print(2, "Vpmgp_size: Invalid mgdisc value (%d)!\n", thee->mgdisc);
236  VASSERT(0);
237  }
238 
239  /* Galerkin or standard coarsening? */
240  switch (thee->mgcoar) { /* NAB TO-DO: This needs to be changed into an enumeration */
241  case 0:
242  if (thee->mgdisc != 0) {
243  Vnm_print(2, "Vpmgp_size: Invalid mgcoar value (%d); must be used with mgdisc 0!\n", thee->mgcoar);
244  VASSERT(0);
245  }
246  num_narrc_oper = 4;
247  break;
248  case 1:
249  if (thee->mgdisc != 0) {
250  Vnm_print(2, "Vpmgp_size: Invalid mgcoar value (%d); must be used with mgdisc 0!\n", thee->mgcoar);
251  VASSERT(0);
252  }
253  num_narrc_oper = 14;
254  break;
255  case 2:
256  num_narrc_oper = 14;
257  break;
258  default:
259  Vnm_print(2, "Vpmgp_size: Invalid mgcoar value (%d)!\n", thee->mgcoar);
260  VASSERT(0);
261  }
262 
263  /* LINPACK storage on coarse grid */
264  switch (thee->mgsolv) { /* NAB TO-DO: This needs to be changed into an enumeration */
265  case 0:
266  n_band = 0;
267  break;
268  case 1:
269  if ( ( (thee->mgcoar == 0) || (thee->mgcoar == 1)) && (thee->mgdisc == 0) ) {
270  num_band = 1 + (thee->nxc-2)*(thee->nyc-2);
271  } else {
272  num_band = 1 + (thee->nxc-2)*(thee->nyc-2) + (thee->nxc-2) + 1;
273  }
274  nc_band = (thee->nxc-2)*(thee->nyc-2)*(thee->nzc-2);
275  n_band = nc_band * num_band;
276  break;
277  default:
278  Vnm_print(2, "Vpmgp_size: Invalid mgsolv value (%d)!\n", thee->mgsolv);
279  VASSERT(0);
280  }
281 
282  /* Real storage parameters */
283  thee->n_rpc = 100*(thee->nlev+1);
284 
285  /* Resulting total required for real storage */
286  thee->nrwk = num_narr*thee->narr + (num_nf + num_nf_oper)*thee->nf + (num_narrc + num_narrc_oper)*thee->narrc + n_band + thee->n_rpc;
287 
288  /* Integer storage parameters */
289  thee->n_iz = 50*(thee->nlev+1);
290  thee->n_ipc = 100*(thee->nlev+1);
291  thee->niwk = thee->n_iz + thee->n_ipc;
292 }
293 
294 VPRIVATE int coarsenThis(int nOld) {
295 
296  int nOut;
297 
298  nOut = (nOld - 1) / 2 + 1;
299 
300  if (((nOut-1)*2) != (nOld-1)) {
301  Vnm_print(2, "Vpmgp_makeCoarse: Warning! The grid dimensions you have chosen are not consistent with the nlev you have specified!\n");
302  Vnm_print(2, "Vpmgp_makeCoarse: This calculation will only work if you are running with mg-dummy type.\n");
303  }
304  if (nOut < 1) {
305  Vnm_print(2, "D'oh! You coarsened the grid below zero! How did you do that?\n");
306  VASSERT(0);
307  }
308 
309  return nOut;
310 }
311 
312 VPUBLIC void Vpmgp_makeCoarse(
313  int numLevel,
314  int nxOld,
315  int nyOld,
316  int nzOld,
317  int *nxNew,
318  int *nyNew,
319  int *nzNew
320  )
321 {
322  int nxtmp, nytmp, nztmp, iLevel;
323 
324  for (iLevel=0; iLevel<numLevel; iLevel++) {
325  nxtmp = *nxNew;
326  nytmp = *nyNew;
327  nztmp = *nzNew;
328  *nxNew = coarsenThis(nxtmp);
329  *nyNew = coarsenThis(nytmp);
330  *nzNew = coarsenThis(nztmp);
331  }
332 
333 
334 }
VPUBLIC void Vpmgp_makeCoarse(int numLevel, int nxOld, int nyOld, int nzOld, int *nxNew, int *nyNew, int *nzNew)
Coarsen the grid by the desired number of levels and determine the resulting numbers of grid points...
Definition: vpmgp.c:312
Definition: vhal.h:211
Contains declarations for class Vpmgp.
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition: vhal.h:563
Contains public data members for Vpmgp class/module.
Definition: vpmgp.h:80
VPUBLIC int Vpmgp_ctor2(Vpmgp *thee, MGparm *mgparm)
FORTRAN stub to construct PMG parameter object and initialize to default values.
Definition: vpmgp.c:93