APBS  1.4.1
mypdec.c
1 
50 #include "mypdec.h"
51 
52 VPUBLIC void Vmypdefinitlpbe(int *tnion, double *tcharge, double *tsconc) {
53 
54  int i;
55 
56  nion = *tnion;
57  if (nion > MAXIONS) {
58  Vnm_print(2, "Vmypde: Warning: Ignoring extra ion species\n");
59  nion = MAXIONS;
60  }
61 
62  for (i=1; i<=nion; i++) {
63  VAT(charge, i) = VAT(tcharge, i);
64  VAT(sconc, i) = VAT(tsconc, i);
65  }
66 }
67 
68 
69 
70 VPUBLIC void Vmypdefinitnpbe(int *tnion, double *tcharge, double *tsconc) {
71 
72  int i;
73 
74  nion = *tnion;
75  if (nion > MAXIONS) {
76  Vnm_print(2, "Vmypde: Warning: Ignoring extra ion species\n");
77  nion = MAXIONS;
78  }
79 
80  for (i=1; i<=nion; i++) {
81  VAT(charge, i) = VAT(tcharge, i);
82  VAT( sconc, i) = VAT( tsconc, i);
83  }
84 }
85 
86 
87 
88 VPUBLIC void Vmypdefinitsmpbe(int *tnion, double *tcharge, double *tsconc,
89  double *smvolume, double *smsize) {
90 
91  int i;
92 
93  WARN_UNTESTED;
94 
95  VABORT_MSG0("Not tested");
96 
97  if (*tnion > 3) {
98  Vnm_print(2, "SMPBE: modified theory handles only three ion species.\n");
99  Vnm_print(2, " Ignoring the rest of the ions!\n");
100  Vnm_print(2, " (mypde.f::mypdefinit)\n");
101  }
102 
103  v1 = VAT(tcharge, 1);
104  v2 = VAT(tcharge, 2);
105  v3 = VAT(tcharge, 3);
106  conc1 = VAT(tsconc, 1);
107  conc2 = VAT(tsconc, 2);
108  conc3 = VAT(tsconc, 3);
109 
110  vol = *smvolume;
111  relSize = *smsize;
112 }
113 
114 
115 
116 VPUBLIC void Vc_vec(double *coef, double *uin, double *uout,
117  int *nx, int *ny, int *nz, int *ipkey) {
118 
119  if (*ipkey == -2) {
120  Vc_vecsmpbe(coef, uin, uout, nx, ny, nz, ipkey);
121  } else {
122  Vc_vecpmg(coef, uin, uout, nx, ny, nz, ipkey);
123  }
124 }
125 
126 
127 
128 VPUBLIC void Vc_vecpmg(double *coef, double *uin, double *uout,
129  int *nx, int *ny, int *nz, int *ipkey) {
130 
131  double zcf2;
132  double zu2;
133  double am_zero;
134  double am_neg;
135  double am_pos;
136  double argument;
137  int ichopped;
138  int ichopped_neg;
139  int ichopped_pos;
140  int iion;
141 
142  int n, i;
143 
144  n = *nx * *ny * *nz;
145 
146  for (i=1; i<=n; i++) {
147  VAT(uout, i) = 0;
148  }
149 
150  for (iion=1; iion<=nion; iion++) {
151 
152  // Assemble the ion-specific coefficient
153  zcf2 = -1.0 * VAT(sconc, iion) * VAT(charge, iion);
154 
155  // Assemble the ion-specific potential value
156  zu2 = -1.0 * VAT(charge, iion);
157 
158  if (*ipkey == 0) {
159  ichopped = 0;
160 
161  #pragma omp parallel for \
162  default(shared) \
163  private(i, ichopped_neg, ichopped_pos, am_zero, am_neg, am_pos, argument) \
164  reduction(+ : ichopped)
165  for (i=1; i<=n; i++) {
166 
167  // am_zero is 0 if coef zero, and 1 if coef nonzero
168  am_zero = VMIN2(ZSMALL, VABS(zcf2 * VAT(coef, i))) * ZLARGE;
169 
170  // am_neg is chopped u if u negative, 0 if u positive
171  am_neg = VMAX2(VMIN2(zu2 * VAT(uin, i), 0.0), SINH_MIN);
172 
173  // am_neg is chopped u if u positive, 0 if u negative
174  am_pos = VMIN2(VMAX2(zu2 * VAT(uin, i), 0.0), SINH_MAX);
175 
176  // Finally determine the function value
177  argument = am_zero * (am_neg + am_pos);
178 
179  VAT(uout, i) = VAT(uout, i) + zcf2 * VAT(coef, i) * exp(argument);
180 
181  // Count chopped values
182  ichopped_neg = (int)(am_neg / SINH_MIN);
183  ichopped_pos = (int)(am_pos / SINH_MAX);
184  ichopped += (int)(floor(am_zero+0.5)) * (ichopped_neg + ichopped_pos);
185  }
186 
187  // Info
188  if (ichopped > 0) {
189  Vnm_print(2, "Vc_vecpmg: trapped exp overflows: %d\n", ichopped);
190  }
191 
192  } else if (*ipkey > 1 && *ipkey % 2 == 1 && *ipkey <= MAXPOLY) {
193 
194  // Polynomial requested
195  Vnm_print(2, "Vc_vecpmg: POLYNOMIAL APPROXIMATION UNAVAILABLE\n");
196  abort();
197  } else {
198 
199  // Return linear approximation !***
200  Vnm_print(2, "Vc_vecpmg: LINEAR APPROXIMATION UNAVAILABLE\n");
201  abort();
202  }
203  }
204 }
205 
206 
207 
208 VPUBLIC void Vc_vecsmpbe(double *coef, double *uin, double *uout,
209  int *nx, int *ny, int *nz, int *ipkey) {
210 
211  int ideg;
212  double zcf2, zu2;
213  double am_zero, am_neg, am_pos;
214  double argument, poly, fact;
215 
216  int ichopped, ichopped_neg, ichopped_pos;
217  int iion;
218  int n, i, ii, ipara, ivect;
219 
220  int nproc = 1;
221 
222  // Added by DG SMPBE variables and common blocks
223  double fracOccA, fracOccB, fracOccC, phi, ionStr;
224  double z1, z2, z3, ca, cb, cc, a, k;
225  double a1_neg, a1_pos, a2_neg, a2_pos;
226  double a3_neg, a3_pos, a1, a2, a3;
227  double f, g, gpark, alpha;
228 
229  WARN_UNTESTED;
230 
231  Vnm_print(2, "Vc_vecsmpbe: v1 = %f\n", v1);
232  Vnm_print(2, "Vc_vecsmpbe: v2 = %f\n", v2);
233  Vnm_print(2, "Vc_vecsmpbe: v3 = %f\n", v3);
234  Vnm_print(2, "Vc_vecsmpbe: conc1 = %f\n", conc1);
235  Vnm_print(2, "Vc_vecsmpbe: conc2 = %f\n", conc2);
236  Vnm_print(2, "Vc_vecsmpbe: conc3 = %f\n", conc3);
237  Vnm_print(2, "Vc_vecsmpbe: vol = %f\n", vol);
238  Vnm_print(2, "Vc_vecsmpbe: relSize = %f\n", relSize);
239 
240  Vnm_print(2, "Vc_vecsmpbe: nion = %d\n", nion);
241 
242  Vnm_print(2, "Vc_vecsmpbe: charge = [");
243  for (i=1; i<=nion; i++)
244  Vnm_print(2, "%f ", VAT(charge, i));
245  Vnm_print(2, "]\n");
246 
247  Vnm_print(2, "Vc_vecsmpbe: sconc = [");
248  for (i=1; i<=nion; i++)
249  Vnm_print(2, "%f ", VAT(sconc, i));
250  Vnm_print(2, "]\n");
251 
252 
253 
254  // Find parallel loops (ipara), remainder (ivect)
255  n = *nx * *ny * *nz;
256  ipara = n / nproc;
257  ivect = n % nproc;
258 
259  for (i=1; i<=n; i++)
260  VAT(uout, i) = 0;
261 
262  // Initialize the chopped counter
263  ichopped = 0;
264 
265  z1 = v1;
266  z2 = v2;
267  z3 = v3;
268  ca = conc1;
269  cb = conc2;
270  cc = conc3;
271  a = vol;
272  k = relSize;
273 
274  if (k - 1 < ZSMALL)
275  Vnm_print(2, "Vc_vecsmpbe: k=1, using special routine\n");
276 
277  // Derived quantities
278  fracOccA = Na * ca * VPOW(a, 3.0);
279  fracOccB = Na * cb * VPOW(a, 3.0);
280  fracOccC = Na * cc * VPOW(a, 3.0);
281 
282  phi = (fracOccA / k) + fracOccB + fracOccC;
283  alpha = (fracOccA / k) / (1 - phi);
284  ionStr = 0.5 * (ca * VPOW(z1, 2.0) + cb * VPOW(z2, 2.0) + cc * VPOW(z3, 2));
285 
286  for (i=1; i<=n; i++) {
287 
288  am_zero = VMIN2(ZSMALL, VABS(VAT(coef, i))) * ZLARGE;
289 
290  // Compute the arguments for exp(-z*u) term
291  a1_neg = VMAX2(VMIN2(-1.0 * z1 * VAT(uin, i), 0.0), SINH_MIN);
292  a1_pos = VMIN2(VMAX2(-1.0 * z1 * VAT(uin, i), 0.0), SINH_MAX);
293 
294  // Compute the arguments for exp(-u) term
295  a2_neg = VMAX2(VMIN2(-1.0 * z2 * VAT(uin, i), 0.0), SINH_MIN);
296  a2_pos = VMIN2(VMAX2(-1.0 * z2 * VAT(uin, i), 0.0), SINH_MAX);
297 
298  // Compute the arguments for exp(u) term
299  a3_neg = VMAX2(VMIN2(-1.0 * z3 * VAT(uin, i), 0.0), SINH_MIN);
300  a3_pos = VMIN2(VMAX2(-1.0 * z3 * VAT(uin, i), 0.0), SINH_MAX);
301 
302  a1 = am_zero * (a1_neg + a1_pos);
303  a2 = am_zero * (a2_neg + a2_pos);
304  a3 = am_zero * (a3_neg + a3_pos);
305 
306  gpark = (1 + alpha * exp(a1)) / (1 + alpha);
307 
308  if (k - 1 < ZSMALL) {
309  f = z1 * ca * exp(a1) + z2 * cb * exp(a2) + z3 * cc * exp(a3);
310  g = 1 - phi + fracOccA * exp(a1)
311  + fracOccB * exp(a2)
312  + fracOccC * exp(a3);
313  } else {
314  f = z1 * ca * exp(a1) * VPOW(gpark, k-1)
315  + z2 * cb * exp(a2)
316  + z3 * cc * exp(a3);
317  g = (1 - phi + fracOccA / k) * VPOW(gpark, k)
318  + fracOccB * exp(a2)
319  + fracOccC * exp(a3);
320  }
321 
322  VAT(uout, i) = -1.0 * VAT(coef, i) * (0.5 / ionStr) * (f / g);
323 
324  // Count chopped values
325  ichopped_neg = (int)((a1_neg + a2_neg+a3_neg) / SINH_MIN);
326  ichopped_pos = (int)((a1_pos + a2_pos+a3_pos) / SINH_MAX);
327  ichopped += (int)floor(am_zero+0.5) * (ichopped_neg + ichopped_pos);
328  }
329 
330  // Info
331  if (ichopped > 0)
332  Vnm_print(2, "Vc_vecsmpbe: trapped exp overflows: %d\n", ichopped);
333 
334 }
335 
336 VPUBLIC void Vdc_vec(double *coef, double *uin, double *uout,
337  int *nx, int *ny, int *nz, int *ipkey) {
338 
339  int i;
340  int n = *nx * *ny * *nz;
341 
342  if(*ipkey == -2) {
343  Vdc_vecsmpbe(coef, uin, uout, nx, ny, nz, ipkey);
344  } else {
345  Vdc_vecpmg(coef, uin, uout, nx, ny, nz, ipkey);
346  }
347 }
348 
349 VPUBLIC void Vdc_vecpmg(double *coef, double *uin, double *uout,
350  int *nx, int *ny, int *nz, int *ipkey) {
351 
352  int ideg, iion;
353  double zcf2, zu2;
354  double am_zero, am_neg, am_pos;
355  double argument, poly, fact;
356 
357  int ichopped, ichopped_neg, ichopped_pos;
358  int n, i;
359 
360  // Find parallel loops (ipara), remainder (ivect)
361  n = *nx * *ny * *nz;
362 
363  for (i=1; i<=n; i++) {
364  VAT(uout, i) = 0.0;
365  }
366 
367  for (iion=1; iion<=nion; iion++) {
368 
369  zcf2 = VAT(sconc, iion) * VAT(charge, iion) * VAT(charge, iion);
370  zu2 = -1.0 * VAT(charge, iion);
371 
372  // Check if full exp requested
373  if (*ipkey == 0) {
374 
375  // Initialize chopped counter
376  ichopped = 0;
377 
378  #pragma omp parallel for \
379  default(shared) \
380  private(i, ichopped_neg, ichopped_pos, \
381  am_zero, am_neg, am_pos, argument) \
382  reduction(+:ichopped)
383  for (i=1; i<=n; i++) {
384 
385  // am_zero is 0 if coef zero, and 1 if coef nonzero
386  am_zero = VMIN2(ZSMALL, VABS(zcf2 * VAT(coef, i))) * ZLARGE;
387 
388  // am_neg is chopped u if u negative, 0 if u positive
389  am_neg = VMAX2(VMIN2(zu2 * VAT(uin, i), 0.0), SINH_MIN);
390 
391  // am_neg is chopped u if u positive, 0 if u negative
392  am_pos = VMIN2(VMAX2(zu2 * VAT(uin, i), 0.0), SINH_MAX);
393 
394  // Finally determine the function value
395  argument = am_zero * (am_neg + am_pos);
396  VAT(uout, i) += zcf2 * VAT(coef, i) * exp( argument );
397 
398  // Count chopped values
399  ichopped_neg = (int)(am_neg / SINH_MIN);
400  ichopped_pos = (int)(am_pos / SINH_MAX);
401  ichopped += (int)floor(am_zero+0.5) * (ichopped_neg + ichopped_pos);
402  }
403 
404  // Info
405  if (ichopped > 0)
406  Vnm_print(2, "Vdc_vec: trapped exp overflows: %d\n", ichopped);
407 
408  } else if ((*ipkey) > 1 && (*ipkey) % 2 == 1 && (*ipkey) <= MAXPOLY) {
409  VABORT_MSG0("Vdc_vec: Polynomial approximation unavailable\n");
410  } else {
411  VABORT_MSG0("Vdc_vec: Linear approximation unavailable\n");
412  }
413  }
414 }
415 
416 
417 
418 VPUBLIC void Vdc_vecsmpbe(double *coef, double *uin, double *uout,
419  int *nx, int *ny, int *nz, int *ipkey) {
420 
421  int ideg, iion;
422  double zcf2, zu2;
423  double am_zero, am_neg, am_pos;
424  double argument, poly, fact;
425  int ichopped, ichopped_neg, ichopped_pos;
426 
427  int n, i, ii;
428  int ipara, ivect;
429 
430  int nproc = 1;
431 
432  // Added by DG SMPBE variables and common blocks
433  double fracOccA, fracOccB, fracOccC, phi, ionStr;
434  double z1, z2, z3, ca, cb, cc, a, k;
435  double a1_neg, a1_pos, a2_neg, a2_pos;
436  double a3_neg, a3_pos, a1, a2, a3;
437  double f, g, fprime, gprime, gpark, alpha;
438 
439  WARN_UNTESTED;
440 
441  // Find parallel loops (ipara), remainder (ivect)
442  n = *nx * *ny * *nz;
443  ipara = n / nproc;
444  ivect = n % nproc;
445 
446  for (i=1; i<=n; i++)
447  VAT(uout, i) = 0.0;
448 
449  // Initialize the chopped counter
450  ichopped = 0;
451 
452  z1 = v1;
453  z2 = v2;
454  z3 = v3;
455  ca = conc1;
456  cb = conc2;
457  cc = conc3;
458  a = vol;
459  k = relSize;
460 
461  if (k - 1 < ZSMALL)
462  Vnm_print(2, "Vdc_vecsmpbe: k=1, using special routine\n");
463 
464  // Derived quantities
465  fracOccA = Na * ca * VPOW(a, 3.0);
466  fracOccB = Na * cb * VPOW(a, 3.0);
467  fracOccC = Na * cc * VPOW(a, 3.0);
468  phi = fracOccA / k + fracOccB + fracOccC;
469  alpha = (fracOccA / k) /(1 - phi);
470  ionStr = 0.5*(ca * VPOW(z1, 2) + cb * VPOW(z2, 2) + cc * VPOW(z3, 2));
471 
472  for (i=1; i<=n; i++) {
473 
474  am_zero = VMIN2(ZSMALL, VABS(VAT(coef, i))) * ZLARGE;
475 
476  // Compute the arguments for exp(-z*u) term
477  a1_neg = VMAX2(VMIN2(-1.0 * z1 * VAT(uin, i), 0.0), SINH_MIN);
478  a1_pos = VMIN2(VMAX2(-1.0 * z1 * VAT(uin, i), 0.0), SINH_MAX);
479 
480  // Compute the arguments for exp(-u) term
481  a2_neg = VMAX2(VMIN2(-1.0 * z2 * VAT(uin, i), 0.0), SINH_MIN);
482  a2_pos = VMIN2(VMAX2(-1.0 * z2 * VAT(uin, i), 0.0), SINH_MAX);
483 
484  // Compute the arguments for exp(u) term
485  a3_neg = VMAX2(VMIN2(-1.0 * z3 * VAT(uin, i), 0.0), SINH_MIN);
486  a3_pos = VMIN2(VMAX2(-1.0 * z3 * VAT(uin, i), 0.0), SINH_MAX);
487 
488  a1 = am_zero * (a1_neg + a1_pos);
489  a2 = am_zero * (a2_neg + a2_pos);
490  a3 = am_zero * (a3_neg + a3_pos);
491 
492  gpark = (1 + alpha * exp(a1)) / (1 + alpha);
493 
494  if (k - 1 < ZSMALL) {
495  f = z1 * ca * exp(a1) + z2 * cb * exp(a2) + z3 * cc * exp(a3);
496  g = 1 - phi + fracOccA * exp(a1)
497  + fracOccB * exp(a2)
498  + fracOccC * exp(a3);
499 
500  fprime =
501  - VPOW(z1, 2) * ca * exp(a1)
502  - VPOW(z2, 2) * cb * exp(a2)
503  - VPOW(z3, 2) * cc * exp(a3);
504 
505  gprime =
506  - z1 * fracOccA * exp(a1)
507  - z2 * fracOccB * exp(a2)
508  - z3 * fracOccC * exp(a3);
509  } else {
510  f = z1 * ca * exp(a1) * VPOW(gpark, k - 1)
511  + z2 * cb * exp(a2)
512  + z3 * cc * exp(a3);
513  g = (1 - phi + fracOccA / k) * VPOW(gpark, k)
514  + fracOccB * exp(a2)
515  + fracOccC * exp(a3);
516 
517  fprime =
518  - VPOW(z1, 2) * ca * exp(a1) * VPOW(gpark, k - 2)
519  * (gpark + (k - 1) * (alpha / (1 + alpha)) * exp(a1))
520  - VPOW(z2, 2) * cb * exp(a2)
521  - VPOW(z3, 2) * cc * exp(a3);
522 
523  gprime =
524  - k * z1 * (alpha / (1 + alpha)) * exp(a1)
525  * (1 - phi + fracOccA / k) * VPOW(gpark, k - 1)
526  - z2 * fracOccB * exp(a2)
527  - z3 * fracOccC * exp(a3);
528 
529  }
530 
531  VAT(uout, i) = -1.0 * VAT(coef, i) * (0.5 / ionStr)
532  * (fprime * g - gprime * f) / VPOW(g, 2.0);
533 
534  // Count chopped values
535  ichopped_neg = (int)((a1_neg + a2_neg + a3_neg) / SINH_MIN);
536  ichopped_pos = (int)((a1_pos + a2_pos + a3_pos) / SINH_MAX);
537  ichopped += (int)floor(am_zero+0.5) * (ichopped_neg + ichopped_pos);
538  }
539 
540  // Info
541  if (ichopped > 0)
542  Vnm_print(2, "Vdc_vecsmpbe: trapped exp overflows: %d\n", ichopped);
543 }
#define SINH_MAX
Used to set the max values acceptable for sinh chopping.
Definition: vhal.h:463
VPUBLIC void Vmypdefinitnpbe(int *tnion, double *tcharge, double *tsconc)
Set up the ionic species to be used in later calculations. This must be called before any other of th...
Definition: mypdec.c:70
int ipkey
Definition: vpmgp.h:109
VPUBLIC void Vc_vecpmg(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
Definition: mypdec.c:128
VPUBLIC void Vc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
Definition: mypdec.c:116
VPUBLIC void Vdc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the derivative of the nonlinearity (vector version)
Definition: mypdec.c:336
VPUBLIC void Vc_vecsmpbe(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
Definition: mypdec.c:208
#define MAXIONS
Specifies the PDE definition for PMG to solve.
Definition: mypdec.h:62
int ny
Definition: vpmgp.h:84
int nx
Definition: vpmgp.h:83
VPUBLIC void Vmypdefinitlpbe(int *tnion, double *tcharge, double *tsconc)
Set up the ionic species to be used in later calculations. This must be called before any other of th...
Definition: mypdec.c:52
#define SINH_MIN
Used to set the min values acceptable for sinh chopping.
Definition: vhal.h:457
VPUBLIC void Vmypdefinitsmpbe(int *tnion, double *tcharge, double *tsconc, double *smvolume, double *smsize)
Set up the ionic species to be used in later calculations. This must be called before any other of th...
Definition: mypdec.c:88
int nz
Definition: vpmgp.h:85