56 #define CENTER(a) a[0] 63 double *rhs,
double bc,
double coeff,
64 double ctr,
int nabor);
72 static void fdaddbc (
int nx,
int ny,
int nz,
int *rp,
int *cval,
73 int *diag,
double *aval,
double *rhs,
double h,
78 #define __FUNC__ "MatGenFDCreate" 91 tmp->
px = tmp->
py = 1;
116 tmp->
a = tmp->
b = tmp->
c = 1.0;
117 tmp->
d = tmp->
e = tmp->
f = 0.0;
118 tmp->
g = tmp->
h = 0.0;
127 tmp->
a = -1 * fabs (tmp->
a);
128 tmp->
b = -1 * fabs (tmp->
b);
129 tmp->
c = -1 * fabs (tmp->
c);
133 tmp->
A = tmp->
B = tmp->
C = tmp->
D = tmp->
E 147 #define __FUNC__ "MatGenFD_Destroy" 157 #define __FUNC__ "MatGenFD_Run" 174 bool debug =
false, striped;
195 int npTest = mg->
px * mg->
py;
201 "numbers don't match: np_dh = %i, px*py*pz = %i",
np,
213 m = mg->
m =
m *
m *
m;
221 mg->
hh = 1.0 / (mg->
px * mg->
cc - 1);
225 sprintf (
msgBuf_dh,
"cc (local grid dimension) = %i", mg->
cc);
247 A->rp = (
int *)
MALLOC_DH ((
m + 1) *
sizeof (int));
286 #define __FUNC__ "generateStriped" 293 int beg_row, end_row;
299 int plane, nodeRemainder;
300 int naborx1, naborx2, nabory1, nabory2;
303 bool applyBdry =
true;
305 double bcx1 = mg->
bcX1;
306 double bcx2 = mg->
bcX2;
307 double bcy1 = mg->
bcY1;
308 double bcy2 = mg->
bcY2;
313 printf_dh (
"@@@ using striped partitioning\n");
324 i = mGlobal / mg->
np;
325 beg_row = i * mg->
id;
326 end_row = beg_row + i;
327 if (mg->
id == mg->
np - 1)
331 mg->
hh = 1.0 / (
m - 1);
332 hhalf = 0.5 * mg->
hh;
335 A->m = end_row - beg_row;
336 A->beg_row = beg_row;
346 fprintf (
logFile,
"generateStriped: beg_row= %i; end_row= %i; m= %i\n",
347 beg_row + 1, end_row + 1,
m);
350 for (row = beg_row; row < end_row; ++row)
352 int localRow = row - beg_row;
356 nodeRemainder = row - (k * plane);
357 j = nodeRemainder /
m;
358 i = nodeRemainder %
m;
362 fprintf (
logFile,
"row= %i x= %i y= %i z= %i\n", row + 1, i, j,
376 cval[idx] = row - plane;
384 nabory1 = cval[idx] = row -
m;
391 naborx1 = cval[idx] = row - 1;
402 naborx2 = cval[idx] = row + 1;
409 nabory2 = cval[idx] = row +
m;
418 cval[idx] = row + plane;
429 int offset = rp[localRow - 1];
430 int len = rp[localRow] - rp[localRow - 1];
437 coeff = mg->
A (mg->
a, i + hhalf, j, k);
438 ctr = mg->
A (mg->
a, i - hhalf, j, k);
440 &(rhs[localRow - 1]), bcx1, coeff, ctr,
443 else if (i == nx - 1)
445 coeff = mg->
A (mg->
a, i - hhalf, j, k);
446 ctr = mg->
A (mg->
a, i + hhalf, j, k);
448 &(rhs[localRow - 1]), bcx2, coeff, ctr,
453 coeff = mg->
B (mg->
b, i, j + hhalf, k);
454 ctr = mg->
B (mg->
b, i, j - hhalf, k);
456 &(rhs[localRow - 1]), bcy1, coeff, ctr,
459 else if (j == ny - 1)
461 coeff = mg->
B (mg->
b, i, j - hhalf, k);
462 ctr = mg->
B (mg->
b, i, j + hhalf, k);
464 &(rhs[localRow - 1]), bcy2, coeff, ctr,
482 const int nx,
const int ny,
const int nz,
int P,
int Q)
485 int lowerx, lowery, lowerz;
507 id = r *
P * Q + q *
P + p;
518 startrow =
id * (nx * ny * nz);
527 return startrow + nx * ny * (z - lowerz) + nx * (y - lowery) + (x -
532 return startrow + nx * (y - lowery) + (x - lowerx);
543 double hhalf =
h * 0.5;
552 for (k = 0; k < 8; ++k)
556 coeff =
g->A (
g->a, x + hhalf, y, z);
560 coeff =
g->A (
g->a, x - hhalf, y, z);
564 coeff =
g->D (
g->d, x, y, z) * hhalf;
569 coeff =
g->B (
g->b, x, y + hhalf, z);
573 coeff =
g->B (
g->b, x, y - hhalf, z);
577 coeff =
g->E (
g->e, x, y, z) * hhalf;
584 coeff =
g->C (
g->c, x, y, z + hhalf);
588 coeff =
g->C (
g->c, x, y, z - hhalf);
592 coeff =
g->F (
g->f, x, y, z) * hhalf;
598 coeff =
g->G (
g->g, x, y, z);
606 konstant (
double coeff,
double x,
double y,
double z)
612 e2_xy (
double coeff,
double x,
double y,
double z)
614 return exp (coeff * x * y);
617 double boxThreeD (
double coeff,
double x,
double y,
double z);
625 box_1 (
double coeff,
double x,
double y,
double z)
627 static bool setup =
false;
628 double retval = coeff;
664 if (x > ax1 && x < ax2 && y > ay1 && y < ay2)
666 retval = dd1 * coeff;
670 if (x > bx1 && x < bx2 && y > by1 && y < by2)
672 retval = dd2 * coeff;
676 if (x > cx1 && x < cx2 && y > cy1 && y < cy2)
678 retval = dd3 * coeff;
687 static bool setup =
false;
688 double retval = coeff;
691 static double dd1 = 100;
694 static double x1 = .2, x2 = .8;
695 static double y1 = .3, y2 = .7;
696 static double z1 = .4, z2 = .6;
706 if (x > x1 && x < x2 && y > y1 && y < y2 && z > z1 && z < z2)
708 retval = dd1 * coeff;
716 box_1 (
double coeff,
double x,
double y,
double z)
718 static double x1, x2, y1, y2;
719 static double d1, d2;
746 if (x > x1 && x < x2 && y > y1 && y < y2)
759 box_2 (
double coeff,
double x,
double y,
double z)
762 static double d1, d2;
775 if (x < .5 && y < .5)
777 if (x > .5 && y > .5)
785 #define __FUNC__ "generateBlocked" 797 int nx =
cc, ny =
cc, nz =
cc;
798 int lowerx, upperx, lowery, uppery, lowerz, upperz;
802 int idx = 0, localRow = 0;
803 int naborx1, naborx2, nabory1, nabory2, naborz1, naborz2;
806 double hhalf = 0.5 * mg->
hh;
807 double bcx1 = mg->
bcX1;
808 double bcx2 = mg->
bcX2;
809 double bcy1 = mg->
bcY1;
810 double bcy2 = mg->
bcY2;
825 q = ((
id - p) /
px) %
py;
826 r = (
id - p -
px * q) / (
px *
py);
831 "this proc's position in subdomain grid: p= %i q= %i r= %i",
840 upperx = lowerx + nx;
842 uppery = lowery + ny;
844 upperz = lowerz + nz;
848 sprintf (
msgBuf_dh,
"local grid parameters: lowerx= %i upperx= %i",
851 sprintf (
msgBuf_dh,
"local grid parameters: lowery= %i uppery= %i",
854 sprintf (
msgBuf_dh,
"local grid parameters: lowerz= %i upperz= %i",
859 startRow = mg->
first;
862 for (z = lowerz; z < upperz; z++)
864 for (y = lowery; y < uppery; y++)
866 for (x = lowerx; x < upperx; x++)
871 fprintf (
logFile,
"row= %i x= %i y= %i z= %i\n",
872 localRow + startRow + 1, x, y, z);
914 cval[idx] = localRow + startRow;
960 int globalRow = localRow + startRow - 1;
961 int offset = rp[localRow - 1];
962 int len = rp[localRow] - rp[localRow - 1];
969 coeff = mg->
A (mg->
a, x + hhalf, y, z);
970 ctr = mg->
A (mg->
a, x - hhalf, y, z);
973 &(rhs[localRow - 1]), bcx1, coeff,
976 else if (x == nx *
px - 1)
978 coeff = mg->
A (mg->
a, x - hhalf, y, z);
979 ctr = mg->
A (mg->
a, x + hhalf, y, z);
982 &(rhs[localRow - 1]), bcx2, coeff,
987 coeff = mg->
B (mg->
b, x, y + hhalf, z);
988 ctr = mg->
B (mg->
b, x, y - hhalf, z);
991 &(rhs[localRow - 1]), bcy1, coeff,
994 else if (y == ny *
py - 1)
996 coeff = mg->
B (mg->
b, x, y - hhalf, z);
997 ctr = mg->
B (mg->
b, x, y + hhalf, z);
1000 &(rhs[localRow - 1]), bcy2, coeff,
1007 coeff = mg->
B (mg->
b, x, y, z + hhalf);
1008 ctr = mg->
B (mg->
b, x, y, z - hhalf);
1011 &(rhs[localRow - 1]), bcy1,
1012 coeff, ctr, naborz2);
1014 else if (z == nz * nx - 1)
1016 coeff = mg->
B (mg->
b, x, y, z - hhalf);
1017 ctr = mg->
B (mg->
b, x, y, z + hhalf);
1020 &(rhs[localRow - 1]), bcy1,
1021 coeff, ctr, naborz1);
1032 #define __FUNC__ "setBoundary_private" 1035 double *rhs,
double bc,
double coeff,
double ctr,
1045 for (i = 0; i < len; ++i)
1047 if (cval[i] == node)
1063 for (i = 0; i < len; ++i)
1066 if (cval[i] == node)
1068 aval[i] += (ctr - coeff);
1071 else if (cval[i] == nabor)
1073 aval[i] = 2.0 * coeff;
double konstant(double coeff, double x, double y, double z)
double(* B)(double coeff, double x, double y, double z)
UF_long CHOLMOD() nnz(cholmod_sparse *A, cholmod_common *Common)
double box_1(double coeff, double x, double y, double z)
static void generateStriped(MatGenFD mg, int *rp, int *cval, double *aval, Mat_dh A, Vec_dh b)
static void setBoundary_private(int node, int *cval, double *aval, int len, double *rhs, double bc, double coeff, double ctr, int nabor)
void Vec_dhInit(Vec_dh v, int size)
bool Parser_dhHasSwitch(Parser_dh p, char *s)
double(* F)(double coeff, double x, double y, double z)
void MatGenFD_Run(MatGenFD mg, int id, int np, Mat_dh *AOut, Vec_dh *rhsOut)
double(* C)(double coeff, double x, double y, double z)
double(* A)(double coeff, double x, double y, double z)
void printf_dh(char *fmt,...)
bool Parser_dhReadDouble(Parser_dh p, char *in, double *out)
double(* G)(double coeff, double x, double y, double z)
bool Parser_dhReadInt(Parser_dh p, char *in, int *out)
double(* E)(double coeff, double x, double y, double z)
int rownum(const bool threeD, const int x, const int y, const int z, const int nx, const int ny, const int nz, int P, int Q)
void Mat_dhCreate(Mat_dh *mat)
static void getstencil(MatGenFD g, int ix, int iy, int iz)
double(* D)(double coeff, double x, double y, double z)
void MatGenFD_Create(MatGenFD *mg)
static void generateBlocked(MatGenFD mg, int *rp, int *cval, double *aval, Mat_dh A, Vec_dh b)
double e2_xy(double coeff, double x, double y, double z)
void MatGenFD_Destroy(MatGenFD mg)
double box_2(double coeff, double x, double y, double z)
double(* H)(double coeff, double x, double y, double z)
void Vec_dhCreate(Vec_dh *v)
char msgBuf_dh[MSG_BUF_SIZE_DH]
double boxThreeD(double coeff, double x, double y, double z)