63 Vnm_print(2,
"Vmypde: Warning: Ignoring extra ion species\n");
67 for (i=1; i<=nion; i++) {
68 VAT(charge, i) = VAT(tcharge, i);
69 VAT(sconc, i) = VAT(tsconc, i);
81 Vnm_print(2,
"Vmypde: Warning: Ignoring extra ion species\n");
85 for (i=1; i<=nion; i++) {
86 VAT(charge, i) = VAT(tcharge, i);
87 VAT( sconc, i) = VAT( tsconc, i);
94 double *smvolume,
double *smsize) {
100 VABORT_MSG0(
"Not tested");
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");
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);
121VPUBLIC
void Vc_vec(
double *coef,
double *uin,
double *uout,
122 int *nx,
int *ny,
int *nz,
int *ipkey) {
127 Vc_vecpmg(coef, uin, uout, nx, ny, nz, ipkey);
133VPUBLIC
void Vc_vecpmg(
double *coef,
double *uin,
double *uout,
134 int *nx,
int *ny,
int *nz,
int *ipkey) {
151 for (i=1; i<=n; i++) {
155 for (iion=1; iion<=nion; iion++) {
158 zcf2 = -1.0 * VAT(sconc, iion) * VAT(charge, iion);
161 zu2 = -1.0 * VAT(charge, iion);
166 #pragma omp parallel for \
168 private(i, ichopped_neg, ichopped_pos, am_zero, am_neg, am_pos, argument) \
169 reduction(+ : ichopped)
170 for (i=1; i<=n; i++) {
173 am_zero = VMIN2(ZSMALL, VABS(zcf2 * VAT(coef, i))) * ZLARGE;
176 am_neg = VMAX2(VMIN2(zu2 * VAT(uin, i), 0.0),
SINH_MIN);
179 am_pos = VMIN2(VMAX2(zu2 * VAT(uin, i), 0.0),
SINH_MAX);
182 argument = am_zero * (am_neg + am_pos);
184 VAT(uout, i) = VAT(uout, i) + zcf2 * VAT(coef, i) * exp(argument);
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);
194 Vnm_print(2,
"Vc_vecpmg: trapped exp overflows: %d\n", ichopped);
197 }
else if (*ipkey > 1 && *ipkey % 2 == 1 && *ipkey <= MAXPOLY) {
200 Vnm_print(2,
"Vc_vecpmg: POLYNOMIAL APPROXIMATION UNAVAILABLE\n");
205 Vnm_print(2,
"Vc_vecpmg: LINEAR APPROXIMATION UNAVAILABLE\n");
214 int *nx,
int *ny,
int *nz,
int *ipkey) {
218 double am_zero, am_neg, am_pos;
219 double argument, poly, fact;
221 int ichopped, ichopped_neg, ichopped_pos;
223 int n, i, ii, ipara, ivect;
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;
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);
245 Vnm_print(2,
"Vc_vecsmpbe: nion = %d\n", nion);
247 Vnm_print(2,
"Vc_vecsmpbe: charge = [");
248 for (i=1; i<=nion; i++)
249 Vnm_print(2,
"%f ", VAT(charge, i));
252 Vnm_print(2,
"Vc_vecsmpbe: sconc = [");
253 for (i=1; i<=nion; i++)
254 Vnm_print(2,
"%f ", VAT(sconc, i));
280 Vnm_print(2,
"Vc_vecsmpbe: k=1, using special routine\n");
283 fracOccA = Na * ca * VPOW(a, 3.0);
284 fracOccB = Na * cb * VPOW(a, 3.0);
285 fracOccC = Na * cc * VPOW(a, 3.0);
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));
291 for (i=1; i<=n; i++) {
293 am_zero = VMIN2(ZSMALL, VABS(VAT(coef, i))) * ZLARGE;
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);
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);
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);
307 a1 = am_zero * (a1_neg + a1_pos);
308 a2 = am_zero * (a2_neg + a2_pos);
309 a3 = am_zero * (a3_neg + a3_pos);
311 gpark = (1 + alpha * exp(a1)) / (1 + alpha);
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)
317 + fracOccC * exp(a3);
319 f = z1 * ca * exp(a1) * VPOW(gpark, k-1)
322 g = (1 - phi + fracOccA / k) * VPOW(gpark, k)
324 + fracOccC * exp(a3);
327 VAT(uout, i) = -1.0 * VAT(coef, i) * (0.5 / ionStr) * (f / g);
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);
337 Vnm_print(2,
"Vc_vecsmpbe: trapped exp overflows: %d\n", ichopped);
341VPUBLIC
void Vdc_vec(
double *coef,
double *uin,
double *uout,
342 int *nx,
int *ny,
int *nz,
int *ipkey) {
345 int n = *nx * *ny * *nz;
348 Vdc_vecsmpbe(coef, uin, uout, nx, ny, nz, ipkey);
350 Vdc_vecpmg(coef, uin, uout, nx, ny, nz, ipkey);
354VPUBLIC
void Vdc_vecpmg(
double *coef,
double *uin,
double *uout,
355 int *nx,
int *ny,
int *nz,
int *ipkey) {
359 double am_zero, am_neg, am_pos;
360 double argument, poly, fact;
362 int ichopped, ichopped_neg, ichopped_pos;
368 for (i=1; i<=n; i++) {
372 for (iion=1; iion<=nion; iion++) {
374 zcf2 = VAT(sconc, iion) * VAT(charge, iion) * VAT(charge, iion);
375 zu2 = -1.0 * VAT(charge, iion);
383 #pragma omp parallel for \
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++) {
391 am_zero = VMIN2(ZSMALL, VABS(zcf2 * VAT(coef, i))) * ZLARGE;
394 am_neg = VMAX2(VMIN2(zu2 * VAT(uin, i), 0.0),
SINH_MIN);
397 am_pos = VMIN2(VMAX2(zu2 * VAT(uin, i), 0.0),
SINH_MAX);
400 argument = am_zero * (am_neg + am_pos);
401 VAT(uout, i) += zcf2 * VAT(coef, i) * exp( argument );
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);
411 Vnm_print(2,
"Vdc_vec: trapped exp overflows: %d\n", ichopped);
413 }
else if ((*ipkey) > 1 && (*ipkey) % 2 == 1 && (*ipkey) <= MAXPOLY) {
414 VABORT_MSG0(
"Vdc_vec: Polynomial approximation unavailable\n");
416 VABORT_MSG0(
"Vdc_vec: Linear approximation unavailable\n");
423VPUBLIC
void Vdc_vecsmpbe(
double *coef,
double *uin,
double *uout,
424 int *nx,
int *ny,
int *nz,
int *ipkey) {
428 double am_zero, am_neg, am_pos;
429 double argument, poly, fact;
430 int ichopped, ichopped_neg, ichopped_pos;
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;
467 Vnm_print(2,
"Vdc_vecsmpbe: k=1, using special routine\n");
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));
477 for (i=1; i<=n; i++) {
479 am_zero = VMIN2(ZSMALL, VABS(VAT(coef, i))) * ZLARGE;
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);
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);
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);
493 a1 = am_zero * (a1_neg + a1_pos);
494 a2 = am_zero * (a2_neg + a2_pos);
495 a3 = am_zero * (a3_neg + a3_pos);
497 gpark = (1 + alpha * exp(a1)) / (1 + alpha);
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)
503 + fracOccC * exp(a3);
506 - VPOW(z1, 2) * ca * exp(a1)
507 - VPOW(z2, 2) * cb * exp(a2)
508 - VPOW(z3, 2) * cc * exp(a3);
511 - z1 * fracOccA * exp(a1)
512 - z2 * fracOccB * exp(a2)
513 - z3 * fracOccC * exp(a3);
515 f = z1 * ca * exp(a1) * VPOW(gpark, k - 1)
518 g = (1 - phi + fracOccA / k) * VPOW(gpark, k)
520 + fracOccC * exp(a3);
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);
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);
536 VAT(uout, i) = -1.0 * VAT(coef, i) * (0.5 / ionStr)
537 * (fprime * g - gprime * f) / VPOW(g, 2.0);
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);
547 Vnm_print(2,
"Vdc_vecsmpbe: trapped exp overflows: %d\n", ichopped);
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 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...
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 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_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
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 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...
#define SINH_MAX
Used to set the max values acceptable for sinh chopping.
#define SINH_MIN
Used to set the min values acceptable for sinh chopping.