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