58 Vnm_print(2,
"Vmypde: Warning: Ignoring extra ion species\n");
62 for (i=1; i<=nion; i++) {
63 VAT(charge, i) = VAT(tcharge, i);
64 VAT(sconc, i) = VAT(tsconc, i);
76 Vnm_print(2,
"Vmypde: Warning: Ignoring extra ion species\n");
80 for (i=1; i<=nion; i++) {
81 VAT(charge, i) = VAT(tcharge, i);
82 VAT( sconc, i) = VAT( tsconc, i);
89 double *smvolume,
double *smsize) {
95 VABORT_MSG0(
"Not tested");
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");
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);
116 VPUBLIC
void Vc_vec(
double *coef,
double *uin,
double *uout,
128 VPUBLIC
void Vc_vecpmg(
double *coef,
double *uin,
double *uout,
146 for (i=1; i<=n; i++) {
150 for (iion=1; iion<=nion; iion++) {
153 zcf2 = -1.0 * VAT(sconc, iion) * VAT(charge, iion);
156 zu2 = -1.0 * VAT(charge, iion);
161 #pragma omp parallel for \ 163 private(i, ichopped_neg, ichopped_pos, am_zero, am_neg, am_pos, argument) \ 164 reduction(+ : ichopped) 165 for (i=1; i<=n; i++) {
168 am_zero = VMIN2(ZSMALL, VABS(zcf2 * VAT(coef, i))) * ZLARGE;
171 am_neg = VMAX2(VMIN2(zu2 * VAT(uin, i), 0.0),
SINH_MIN);
174 am_pos = VMIN2(VMAX2(zu2 * VAT(uin, i), 0.0),
SINH_MAX);
177 argument = am_zero * (am_neg + am_pos);
179 VAT(uout, i) = VAT(uout, i) + zcf2 * VAT(coef, i) * exp(argument);
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);
189 Vnm_print(2,
"Vc_vecpmg: trapped exp overflows: %d\n", ichopped);
195 Vnm_print(2,
"Vc_vecpmg: POLYNOMIAL APPROXIMATION UNAVAILABLE\n");
200 Vnm_print(2,
"Vc_vecpmg: LINEAR APPROXIMATION UNAVAILABLE\n");
213 double am_zero, am_neg, am_pos;
214 double argument, poly, fact;
216 int ichopped, ichopped_neg, ichopped_pos;
218 int n, i, ii, ipara, ivect;
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;
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);
240 Vnm_print(2,
"Vc_vecsmpbe: nion = %d\n", nion);
242 Vnm_print(2,
"Vc_vecsmpbe: charge = [");
243 for (i=1; i<=nion; i++)
244 Vnm_print(2,
"%f ", VAT(charge, i));
247 Vnm_print(2,
"Vc_vecsmpbe: sconc = [");
248 for (i=1; i<=nion; i++)
249 Vnm_print(2,
"%f ", VAT(sconc, i));
275 Vnm_print(2,
"Vc_vecsmpbe: k=1, using special routine\n");
278 fracOccA = Na * ca * VPOW(a, 3.0);
279 fracOccB = Na * cb * VPOW(a, 3.0);
280 fracOccC = Na * cc * VPOW(a, 3.0);
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));
286 for (i=1; i<=n; i++) {
288 am_zero = VMIN2(ZSMALL, VABS(VAT(coef, i))) * ZLARGE;
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);
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);
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);
302 a1 = am_zero * (a1_neg + a1_pos);
303 a2 = am_zero * (a2_neg + a2_pos);
304 a3 = am_zero * (a3_neg + a3_pos);
306 gpark = (1 + alpha * exp(a1)) / (1 + alpha);
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)
312 + fracOccC * exp(a3);
314 f = z1 * ca * exp(a1) * VPOW(gpark, k-1)
317 g = (1 - phi + fracOccA / k) * VPOW(gpark, k)
319 + fracOccC * exp(a3);
322 VAT(uout, i) = -1.0 * VAT(coef, i) * (0.5 / ionStr) * (f / g);
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);
332 Vnm_print(2,
"Vc_vecsmpbe: trapped exp overflows: %d\n", ichopped);
336 VPUBLIC
void Vdc_vec(
double *coef,
double *uin,
double *uout,
349 VPUBLIC
void Vdc_vecpmg(
double *coef,
double *uin,
double *uout,
354 double am_zero, am_neg, am_pos;
355 double argument, poly, fact;
357 int ichopped, ichopped_neg, ichopped_pos;
363 for (i=1; i<=n; i++) {
367 for (iion=1; iion<=nion; iion++) {
369 zcf2 = VAT(sconc, iion) * VAT(charge, iion) * VAT(charge, iion);
370 zu2 = -1.0 * VAT(charge, iion);
378 #pragma omp parallel for \ 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++) {
386 am_zero = VMIN2(ZSMALL, VABS(zcf2 * VAT(coef, i))) * ZLARGE;
389 am_neg = VMAX2(VMIN2(zu2 * VAT(uin, i), 0.0),
SINH_MIN);
392 am_pos = VMIN2(VMAX2(zu2 * VAT(uin, i), 0.0),
SINH_MAX);
395 argument = am_zero * (am_neg + am_pos);
396 VAT(uout, i) += zcf2 * VAT(coef, i) * exp( argument );
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);
406 Vnm_print(2,
"Vdc_vec: trapped exp overflows: %d\n", ichopped);
409 VABORT_MSG0(
"Vdc_vec: Polynomial approximation unavailable\n");
411 VABORT_MSG0(
"Vdc_vec: Linear approximation unavailable\n");
418 VPUBLIC
void Vdc_vecsmpbe(
double *coef,
double *uin,
double *uout,
423 double am_zero, am_neg, am_pos;
424 double argument, poly, fact;
425 int ichopped, ichopped_neg, ichopped_pos;
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;
462 Vnm_print(2,
"Vdc_vecsmpbe: k=1, using special routine\n");
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));
472 for (i=1; i<=n; i++) {
474 am_zero = VMIN2(ZSMALL, VABS(VAT(coef, i))) * ZLARGE;
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);
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);
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);
488 a1 = am_zero * (a1_neg + a1_pos);
489 a2 = am_zero * (a2_neg + a2_pos);
490 a3 = am_zero * (a3_neg + a3_pos);
492 gpark = (1 + alpha * exp(a1)) / (1 + alpha);
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)
498 + fracOccC * exp(a3);
501 - VPOW(z1, 2) * ca * exp(a1)
502 - VPOW(z2, 2) * cb * exp(a2)
503 - VPOW(z3, 2) * cc * exp(a3);
506 - z1 * fracOccA * exp(a1)
507 - z2 * fracOccB * exp(a2)
508 - z3 * fracOccC * exp(a3);
510 f = z1 * ca * exp(a1) * VPOW(gpark, k - 1)
513 g = (1 - phi + fracOccA / k) * VPOW(gpark, k)
515 + fracOccC * exp(a3);
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);
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);
531 VAT(uout, i) = -1.0 * VAT(coef, i) * (0.5 / ionStr)
532 * (fprime * g - gprime * f) / VPOW(g, 2.0);
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);
542 Vnm_print(2,
"Vdc_vecsmpbe: trapped exp overflows: %d\n", ichopped);
#define SINH_MAX
Used to set the max values acceptable for sinh chopping.
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...
VPUBLIC void Vc_vecpmg(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
VPUBLIC void Vc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
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)
VPUBLIC void Vc_vecsmpbe(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
#define MAXIONS
Specifies the PDE definition for PMG to solve.
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...
#define SINH_MIN
Used to set the min values acceptable for sinh chopping.
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...