APBS 3.0.0
Loading...
Searching...
No Matches
mypdec.c
1
55#include "mypdec.h"
56
57VPUBLIC 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
75VPUBLIC 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
93VPUBLIC 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
121VPUBLIC 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
133VPUBLIC 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
213VPUBLIC 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
341VPUBLIC 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
354VPUBLIC 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
423VPUBLIC 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