IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
MatGenFD.c
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "MatGenFD.h"
44 #include "Mat_dh.h"
45 #include "Vec_dh.h"
46 #include "Parser_dh.h"
47 #include "Mem_dh.h"
48 /* #include "graphColor_dh.h" */
49 
50 static bool isThreeD;
51 
52  /* handles for values in the 5-point (2D) or 7-point (for 3D) stencil */
53 #define FRONT(a) a[5]
54 #define SOUTH(a) a[3]
55 #define WEST(a) a[1]
56 #define CENTER(a) a[0]
57 #define EAST(a) a[2]
58 #define NORTH(a) a[4]
59 #define BACK(a) a[6]
60 #define RHS(a) a[7]
61 
62 static void setBoundary_private (int node, int *cval, double *aval, int len,
63  double *rhs, double bc, double coeff,
64  double ctr, int nabor);
65 static void generateStriped (MatGenFD mg, int *rp, int *cval, double *aval,
66  Mat_dh A, Vec_dh b);
67 static void generateBlocked (MatGenFD mg, int *rp, int *cval, double *aval,
68  Mat_dh A, Vec_dh b);
69 static void getstencil (MatGenFD g, int ix, int iy, int iz);
70 
71 #if 0
72 static void fdaddbc (int nx, int ny, int nz, int *rp, int *cval,
73  int *diag, double *aval, double *rhs, double h,
74  MatGenFD mg);
75 #endif
76 
77 #undef __FUNC__
78 #define __FUNC__ "MatGenFDCreate"
79 void
80 MatGenFD_Create (MatGenFD * mg)
81 {
82  START_FUNC_DH
83  struct _matgenfd *tmp =
84  (struct _matgenfd *) MALLOC_DH (sizeof (struct _matgenfd));
85  CHECK_V_ERROR;
86  *mg = tmp;
87 
88  tmp->debug = Parser_dhHasSwitch (parser_dh, "-debug_matgen");
89 
90  tmp->m = 9;
91  tmp->px = tmp->py = 1;
92  tmp->pz = 0;
93  Parser_dhReadInt (parser_dh, "-m", &tmp->m);
94  Parser_dhReadInt (parser_dh, "-px", &tmp->px);
95  Parser_dhReadInt (parser_dh, "-py", &tmp->py);
96  Parser_dhReadInt (parser_dh, "-pz", &tmp->pz);
97 
98  if (tmp->px < 1)
99  tmp->px = 1;
100  if (tmp->py < 1)
101  tmp->py = 1;
102  if (tmp->pz < 0)
103  tmp->pz = 0;
104  tmp->threeD = false;
105  if (tmp->pz)
106  {
107  tmp->threeD = true;
108  }
109  else
110  {
111  tmp->pz = 1;
112  }
113  if (Parser_dhHasSwitch (parser_dh, "-threeD"))
114  tmp->threeD = true;
115 
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;
119 
120  Parser_dhReadDouble (parser_dh, "-dx", &tmp->a);
121  Parser_dhReadDouble (parser_dh, "-dy", &tmp->b);
122  Parser_dhReadDouble (parser_dh, "-dz", &tmp->c);
123  Parser_dhReadDouble (parser_dh, "-cx", &tmp->d);
124  Parser_dhReadDouble (parser_dh, "-cy", &tmp->e);
125  Parser_dhReadDouble (parser_dh, "-cz", &tmp->f);
126 
127  tmp->a = -1 * fabs (tmp->a);
128  tmp->b = -1 * fabs (tmp->b);
129  tmp->c = -1 * fabs (tmp->c);
130 
131  tmp->allocateMem = true;
132 
133  tmp->A = tmp->B = tmp->C = tmp->D = tmp->E
134  = tmp->F = tmp->G = tmp->H = konstant;
135 
136  tmp->bcX1 = tmp->bcX2 = tmp->bcY1 = tmp->bcY2 = tmp->bcZ1 = tmp->bcZ2 = 0.0;
137  Parser_dhReadDouble (parser_dh, "-bcx1", &tmp->bcX1);
138  Parser_dhReadDouble (parser_dh, "-bcx2", &tmp->bcX2);
139  Parser_dhReadDouble (parser_dh, "-bcy1", &tmp->bcY1);
140  Parser_dhReadDouble (parser_dh, "-bcy2", &tmp->bcY2);
141  Parser_dhReadDouble (parser_dh, "-bcz1", &tmp->bcZ1);
142  Parser_dhReadDouble (parser_dh, "-bcz2", &tmp->bcZ2);
143 END_FUNC_DH}
144 
145 
146 #undef __FUNC__
147 #define __FUNC__ "MatGenFD_Destroy"
148 void
149 MatGenFD_Destroy (MatGenFD mg)
150 {
151  START_FUNC_DH FREE_DH (mg);
152  CHECK_V_ERROR;
153 END_FUNC_DH}
154 
155 
156 #undef __FUNC__
157 #define __FUNC__ "MatGenFD_Run"
158 void
159 MatGenFD_Run (MatGenFD mg, int id, int np, Mat_dh * AOut, Vec_dh * rhsOut)
160 {
161 /* What this function does:
162  * 0. creates return objects (A and rhs)
163  * 1. computes "nice to have" values;
164  * 2. allocates storage, if required;
165  * 3. calls generateBlocked() or generateStriped().
166  * 4. initializes variable in A and rhs.
167  */
168 
169  START_FUNC_DH Mat_dh A;
170  Vec_dh rhs;
171  bool threeD = mg->threeD;
172  int nnz;
173  int m = mg->m; /* local unknowns */
174  bool debug = false, striped;
175 
176  if (mg->debug && logFile != NULL)
177  debug = true;
178  striped = Parser_dhHasSwitch (parser_dh, "-striped");
179 
180  /* 0. create objects */
181  Mat_dhCreate (AOut);
182  CHECK_V_ERROR;
183  Vec_dhCreate (rhsOut);
184  CHECK_V_ERROR;
185  A = *AOut;
186  rhs = *rhsOut;
187 
188  /* ensure that processor grid contains the same number of
189  nodes as there are processors.
190  */
191  if (!Parser_dhHasSwitch (parser_dh, "-noChecks"))
192  {
193  if (!striped)
194  {
195  int npTest = mg->px * mg->py;
196  if (threeD)
197  npTest *= mg->pz;
198  if (npTest != np)
199  {
200  sprintf (msgBuf_dh,
201  "numbers don't match: np_dh = %i, px*py*pz = %i", np,
202  npTest);
203  SET_V_ERROR (msgBuf_dh);
204  }
205  }
206  }
207 
208  /* 1. compute "nice to have" values */
209  /* each proc's subgrid dimension */
210  mg->cc = m;
211  if (threeD)
212  {
213  m = mg->m = m * m * m;
214  }
215  else
216  {
217  m = mg->m = m * m;
218  }
219 
220  mg->first = id * m;
221  mg->hh = 1.0 / (mg->px * mg->cc - 1);
222 
223  if (debug)
224  {
225  sprintf (msgBuf_dh, "cc (local grid dimension) = %i", mg->cc);
226  SET_INFO (msgBuf_dh);
227  if (threeD)
228  {
229  sprintf (msgBuf_dh, "threeD = true");
230  }
231  else
232  {
233  sprintf (msgBuf_dh, "threeD = false");
234  }
235  SET_INFO (msgBuf_dh);
236  sprintf (msgBuf_dh, "np= %i id= %i", np, id);
237  SET_INFO (msgBuf_dh);
238  }
239 
240  mg->id = id;
241  mg->np = np;
242  nnz = threeD ? m * 7 : m * 5;
243 
244  /* 2. allocate storage */
245  if (mg->allocateMem)
246  {
247  A->rp = (int *) MALLOC_DH ((m + 1) * sizeof (int));
248  CHECK_V_ERROR;
249  A->rp[0] = 0;
250  A->cval = (int *) MALLOC_DH (nnz * sizeof (int));
251  CHECK_V_ERROR A->aval = (double *) MALLOC_DH (nnz * sizeof (double));
252  CHECK_V_ERROR;
253  /* rhs->vals = (double*)MALLOC_DH(m*sizeof(double)); CHECK_V_ERROR; */
254  }
255 
256  /* 4. initialize variables in A and rhs */
257  rhs->n = m;
258  A->m = m;
259  A->n = m * mg->np;
260  A->beg_row = mg->first;
261 
262  /* 3. generate matrix */
263  isThreeD = threeD; /* yuck! used in box_XX() */
264  if (Parser_dhHasSwitch (parser_dh, "-striped"))
265  {
266  generateStriped (mg, A->rp, A->cval, A->aval, A, rhs);
267  CHECK_V_ERROR;
268  }
269  else
270  {
271  generateBlocked (mg, A->rp, A->cval, A->aval, A, rhs);
272  CHECK_V_ERROR;
273  }
274 
275  /* add in bdry conditions */
276  /* only implemented for 2D mats! */
277  if (!threeD)
278  {
279 /* fdaddbc(nx, ny, nz, rp, cval, diag, aval, rhs, h, mg); */
280  }
281 
282 END_FUNC_DH}
283 
284 
285 #undef __FUNC__
286 #define __FUNC__ "generateStriped"
287 void
288 generateStriped (MatGenFD mg, int *rp, int *cval, double *aval, Mat_dh A,
289  Vec_dh b)
290 {
291  START_FUNC_DH int mGlobal;
292  int m = mg->m;
293  int beg_row, end_row;
294  int i, j, k, row;
295  bool threeD = mg->threeD;
296  int idx = 0;
297  double *stencil = mg->stencil;
298  bool debug = false;
299  int plane, nodeRemainder;
300  int naborx1, naborx2, nabory1, nabory2;
301  double *rhs;
302 
303  bool applyBdry = true;
304  double hhalf;
305  double bcx1 = mg->bcX1;
306  double bcx2 = mg->bcX2;
307  double bcy1 = mg->bcY1;
308  double bcy2 = mg->bcY2;
309  /* double bcz1 = mg->bcZ1; */
310  /* double bcz2 = mg->bcZ2; */
311  int nx, ny;
312 
313  printf_dh ("@@@ using striped partitioning\n");
314 
315  if (mg->debug && logFile != NULL)
316  debug = true;
317 
318  /* recompute values (yuck!) */
319  m = 9;
320  Parser_dhReadInt (parser_dh, "-m", &m); /* global grid dimension */
321  mGlobal = m * m; /* global unkknowns */
322  if (threeD)
323  mGlobal *= m;
324  i = mGlobal / mg->np; /* unknowns per processor */
325  beg_row = i * mg->id; /* global number of 1st local row */
326  end_row = beg_row + i;
327  if (mg->id == mg->np - 1)
328  end_row = mGlobal;
329  nx = ny = m;
330 
331  mg->hh = 1.0 / (m - 1);
332  hhalf = 0.5 * mg->hh;
333 
334  A->n = m * m;
335  A->m = end_row - beg_row;
336  A->beg_row = beg_row;
337 
338  Vec_dhInit (b, A->m);
339  CHECK_V_ERROR;
340  rhs = b->vals;
341 
342  plane = m * m;
343 
344  if (debug)
345  {
346  fprintf (logFile, "generateStriped: beg_row= %i; end_row= %i; m= %i\n",
347  beg_row + 1, end_row + 1, m);
348  }
349 
350  for (row = beg_row; row < end_row; ++row)
351  {
352  int localRow = row - beg_row;
353 
354  /* compute current node's position in grid */
355  k = (row / plane);
356  nodeRemainder = row - (k * plane); /* map row to 1st plane */
357  j = nodeRemainder / m;
358  i = nodeRemainder % m;
359 
360  if (debug)
361  {
362  fprintf (logFile, "row= %i x= %i y= %i z= %i\n", row + 1, i, j,
363  k);
364  }
365 
366  /* compute column values and rhs entry for the current node */
367  getstencil (mg, i, j, k);
368 
369  /* only homogenous Dirichlet boundary conditions presently supported */
370 
371  /* down plane */
372  if (threeD)
373  {
374  if (k > 0)
375  {
376  cval[idx] = row - plane;
377  aval[idx++] = BACK (stencil);
378  }
379  }
380 
381  /* south */
382  if (j > 0)
383  {
384  nabory1 = cval[idx] = row - m;
385  aval[idx++] = SOUTH (stencil);
386  }
387 
388  /* west */
389  if (i > 0)
390  {
391  naborx1 = cval[idx] = row - 1;
392  aval[idx++] = WEST (stencil);
393  }
394 
395  /* center node */
396  cval[idx] = row;
397  aval[idx++] = CENTER (stencil);
398 
399  /* east */
400  if (i < m - 1)
401  {
402  naborx2 = cval[idx] = row + 1;
403  aval[idx++] = EAST (stencil);
404  }
405 
406  /* north */
407  if (j < m - 1)
408  {
409  nabory2 = cval[idx] = row + m;
410  aval[idx++] = NORTH (stencil);
411  }
412 
413  /* up plane */
414  if (threeD)
415  {
416  if (k < m - 1)
417  {
418  cval[idx] = row + plane;
419  aval[idx++] = FRONT (stencil);
420  }
421  }
422  rhs[localRow] = 0.0;
423  ++localRow;
424  rp[localRow] = idx;
425 
426  /* apply boundary conditions; only for 2D! */
427  if (!threeD && applyBdry)
428  {
429  int offset = rp[localRow - 1];
430  int len = rp[localRow] - rp[localRow - 1];
431  double ctr, coeff;
432 
433 /* fprintf(logFile, "globalRow = %i; naborx2 = %i\n", row+1, row); */
434 
435  if (i == 0)
436  { /* if x1 */
437  coeff = mg->A (mg->a, i + hhalf, j, k);
438  ctr = mg->A (mg->a, i - hhalf, j, k);
439  setBoundary_private (row, cval + offset, aval + offset, len,
440  &(rhs[localRow - 1]), bcx1, coeff, ctr,
441  naborx2);
442  }
443  else if (i == nx - 1)
444  { /* if x2 */
445  coeff = mg->A (mg->a, i - hhalf, j, k);
446  ctr = mg->A (mg->a, i + hhalf, j, k);
447  setBoundary_private (row, cval + offset, aval + offset, len,
448  &(rhs[localRow - 1]), bcx2, coeff, ctr,
449  naborx1);
450  }
451  else if (j == 0)
452  { /* if y1 */
453  coeff = mg->B (mg->b, i, j + hhalf, k);
454  ctr = mg->B (mg->b, i, j - hhalf, k);
455  setBoundary_private (row, cval + offset, aval + offset, len,
456  &(rhs[localRow - 1]), bcy1, coeff, ctr,
457  nabory2);
458  }
459  else if (j == ny - 1)
460  { /* if y2 */
461  coeff = mg->B (mg->b, i, j - hhalf, k);
462  ctr = mg->B (mg->b, i, j + hhalf, k);
463  setBoundary_private (row, cval + offset, aval + offset, len,
464  &(rhs[localRow - 1]), bcy2, coeff, ctr,
465  nabory1);
466  }
467  }
468  }
469 END_FUNC_DH}
470 
471 
472 /* zero-based
473  (from Edmond Chow)
474 */
475 /*
476  x,y,z - coordinates of row, wrt naturally ordered grid
477  nz, ny, nz - local grid dimensions, wrt 0
478  P, Q - subdomain grid dimensions in x and y directions
479 */
480 int
481 rownum (const bool threeD, const int x, const int y, const int z,
482  const int nx, const int ny, const int nz, int P, int Q)
483 {
484  int p, q, r;
485  int lowerx, lowery, lowerz;
486  int id, startrow;
487 
488 
489  /* compute x,y,z coordinates of subdomain to which
490  this row belongs.
491  */
492  p = x / nx;
493  q = y / ny;
494  r = z / nz;
495 
496 /*
497 if (myid_dh == 0) printf("nx= %i ny= %i nz= %i\n", nx, ny, nz);
498 if (myid_dh == 0) printf("x= %i y= %i z= %i threeD= %i p= %i q= %i r= %i\n",
499  x,y,z,threeD, p,q,r);
500 */
501 
502  /* compute the subdomain (processor) of the subdomain to which
503  this row belongs.
504  */
505  if (threeD)
506  {
507  id = r * P * Q + q * P + p;
508  }
509  else
510  {
511  id = q * P + p;
512  }
513 
514 /* if (myid_dh == 0) printf(" id= %i\n", id);
515 */
516 
517  /* smallest row in the subdomain */
518  startrow = id * (nx * ny * nz);
519 
520  /* x,y, and z coordinates of local grid of unknowns */
521  lowerx = nx * p;
522  lowery = ny * q;
523  lowerz = nz * r;
524 
525  if (threeD)
526  {
527  return startrow + nx * ny * (z - lowerz) + nx * (y - lowery) + (x -
528  lowerx);
529  }
530  else
531  {
532  return startrow + nx * (y - lowery) + (x - lowerx);
533  }
534 }
535 
536 
537 
538 void
539 getstencil (MatGenFD g, int ix, int iy, int iz)
540 {
541  int k;
542  double h = g->hh;
543  double hhalf = h * 0.5;
544  double x = h * ix;
545  double y = h * iy;
546  double z = h * iz;
547  double cntr = 0.0;
548  double *stencil = g->stencil;
549  double coeff;
550  bool threeD = g->threeD;
551 
552  for (k = 0; k < 8; ++k)
553  stencil[k] = 0.0;
554 
555  /* differentiation wrt x */
556  coeff = g->A (g->a, x + hhalf, y, z);
557  EAST (stencil) += coeff;
558  cntr += coeff;
559 
560  coeff = g->A (g->a, x - hhalf, y, z);
561  WEST (stencil) += coeff;
562  cntr += coeff;
563 
564  coeff = g->D (g->d, x, y, z) * hhalf;
565  EAST (stencil) += coeff;
566  WEST (stencil) -= coeff;
567 
568  /* differentiation wrt y */
569  coeff = g->B (g->b, x, y + hhalf, z);
570  NORTH (stencil) += coeff;
571  cntr += coeff;
572 
573  coeff = g->B (g->b, x, y - hhalf, z);
574  SOUTH (stencil) += coeff;
575  cntr += coeff;
576 
577  coeff = g->E (g->e, x, y, z) * hhalf;
578  NORTH (stencil) += coeff;
579  SOUTH (stencil) -= coeff;
580 
581  /* differentiation wrt z */
582  if (threeD)
583  {
584  coeff = g->C (g->c, x, y, z + hhalf);
585  BACK (stencil) += coeff;
586  cntr += coeff;
587 
588  coeff = g->C (g->c, x, y, z - hhalf);
589  FRONT (stencil) += coeff;
590  cntr += coeff;
591 
592  coeff = g->F (g->f, x, y, z) * hhalf;
593  BACK (stencil) += coeff;
594  FRONT (stencil) -= coeff;
595  }
596 
597  /* contribution from function G: */
598  coeff = g->G (g->g, x, y, z);
599  CENTER (stencil) = h * h * coeff - cntr;
600 
601  RHS (stencil) = h * h * g->H (g->h, x, y, z);
602 }
603 
604 
605 double
606 konstant (double coeff, double x, double y, double z)
607 {
608  return coeff;
609 }
610 
611 double
612 e2_xy (double coeff, double x, double y, double z)
613 {
614  return exp (coeff * x * y);
615 }
616 
617 double boxThreeD (double coeff, double x, double y, double z);
618 
619 /* returns diffusivity constant -bd1 if the point
620  (x,y,z) is inside the box whose upper left and
621  lower right points are (-bx1,-by1), (-bx2,-by2);
622  else, returns diffusivity constant -bd2
623 */
624 double
625 box_1 (double coeff, double x, double y, double z)
626 {
627  static bool setup = false;
628  double retval = coeff;
629 
630  /* dffusivity constants */
631  static double dd1 = BOX1_DD;
632  static double dd2 = BOX2_DD;
633  static double dd3 = BOX3_DD;
634 
635  /* boxes */
636  static double ax1 = BOX1_X1, ay1 = BOX1_Y1;
637  static double ax2 = BOX1_X2, ay2 = BOX1_Y2;
638  static double bx1 = BOX2_X1, by1 = BOX2_Y1;
639  static double bx2 = BOX2_X2, by2 = BOX2_Y2;
640  static double cx1 = BOX3_X1, cy1 = BOX3_Y1;
641  static double cx2 = BOX3_X2, cy2 = BOX3_Y2;
642 
643  if (isThreeD)
644  {
645  return (boxThreeD (coeff, x, y, z));
646  }
647 
648 
649  /* 1st time through, parse for dffusivity constants */
650  if (!setup)
651  {
652  dd1 = 0.1;
653  dd2 = 0.1;
654  dd3 = 10;
655  Parser_dhReadDouble (parser_dh, "-dd1", &dd1);
656  Parser_dhReadDouble (parser_dh, "-dd2", &dd2);
657  Parser_dhReadDouble (parser_dh, "-dd3", &dd3);
658  Parser_dhReadDouble (parser_dh, "-box1x1", &cx1);
659  Parser_dhReadDouble (parser_dh, "-box1x2", &cx2);
660  setup = true;
661  }
662 
663  /* determine if point is inside box a */
664  if (x > ax1 && x < ax2 && y > ay1 && y < ay2)
665  {
666  retval = dd1 * coeff;
667  }
668 
669  /* determine if point is inside box b */
670  if (x > bx1 && x < bx2 && y > by1 && y < by2)
671  {
672  retval = dd2 * coeff;
673  }
674 
675  /* determine if point is inside box c */
676  if (x > cx1 && x < cx2 && y > cy1 && y < cy2)
677  {
678  retval = dd3 * coeff;
679  }
680 
681  return retval;
682 }
683 
684 double
685 boxThreeD (double coeff, double x, double y, double z)
686 {
687  static bool setup = false;
688  double retval = coeff;
689 
690  /* dffusivity constants */
691  static double dd1 = 100;
692 
693  /* boxes */
694  static double x1 = .2, x2 = .8;
695  static double y1 = .3, y2 = .7;
696  static double z1 = .4, z2 = .6;
697 
698  /* 1st time through, parse for diffusivity constants */
699  if (!setup)
700  {
701  Parser_dhReadDouble (parser_dh, "-dd1", &dd1);
702  setup = true;
703  }
704 
705  /* determine if point is inside the box */
706  if (x > x1 && x < x2 && y > y1 && y < y2 && z > z1 && z < z2)
707  {
708  retval = dd1 * coeff;
709  }
710 
711  return retval;
712 }
713 
714 #if 0
715 double
716 box_1 (double coeff, double x, double y, double z)
717 {
718  static double x1, x2, y1, y2;
719  static double d1, d2;
720  bool setup = false;
721  double retval;
722 
723  /* 1st time through, parse for constants and
724  bounding box definition
725  */
726  if (!setup)
727  {
728  x1 = .25;
729  x2 = .75;
730  y1 = .25;
731  y2 = .75;
732  d1 = 1;
733  d2 = 2;
734  Parser_dhReadDouble (parser_dh, "-bx1", &x1);
735  Parser_dhReadDouble (parser_dh, "-bx2", &x2);
736  Parser_dhReadDouble (parser_dh, "-by1", &y1);
737  Parser_dhReadDouble (parser_dh, "-by2", &y2);
738  Parser_dhReadDouble (parser_dh, "-bd1", &d1);
739  Parser_dhReadDouble (parser_dh, "-bd2", &d2);
740  setup = true;
741  }
742 
743  retval = d2;
744 
745  /* determine if point is inside box */
746  if (x > x1 && x < x2 && y > y1 && y < y2)
747  {
748  retval = d1;
749  }
750 
751  return -1 * retval;
752 }
753 #endif
754 
755 /* divide square into 4 quadrants; return one of
756  2 constants depending on the quadrant (checkerboard)
757 */
758 double
759 box_2 (double coeff, double x, double y, double z)
760 {
761  bool setup = false;
762  static double d1, d2;
763  double retval;
764 
765  if (!setup)
766  {
767  d1 = 1;
768  d2 = 2;
769  Parser_dhReadDouble (parser_dh, "-bd1", &d1);
770  Parser_dhReadDouble (parser_dh, "-bd2", &d2);
771  }
772 
773  retval = d2;
774 
775  if (x < .5 && y < .5)
776  retval = d1;
777  if (x > .5 && y > .5)
778  retval = d1;
779 
780  return -1 * retval;
781 }
782 
783 
784 #undef __FUNC__
785 #define __FUNC__ "generateBlocked"
786 void
787 generateBlocked (MatGenFD mg, int *rp, int *cval, double *aval, Mat_dh A,
788  Vec_dh b)
789 {
790  START_FUNC_DH bool applyBdry = true;
791  double *stencil = mg->stencil;
792  int id = mg->id;
793  bool threeD = mg->threeD;
794  int px = mg->px, py = mg->py, pz = mg->pz; /* processor grid dimensions */
795  int p, q, r; /* this proc's position in processor grid */
796  int cc = mg->cc; /* local grid dimension (grid of unknowns) */
797  int nx = cc, ny = cc, nz = cc;
798  int lowerx, upperx, lowery, uppery, lowerz, upperz;
799  int startRow;
800  int x, y, z;
801  bool debug = false;
802  int idx = 0, localRow = 0; /* nabor; */
803  int naborx1, naborx2, nabory1, nabory2, naborz1, naborz2;
804  double *rhs;
805 
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;
811  /* double bcz1 = mg->bcZ1; */
812  /* double bcz2 = mg->bcZ2; */
813 
814  Vec_dhInit (b, A->m);
815  CHECK_V_ERROR;
816  rhs = b->vals;
817 
818  if (mg->debug && logFile != NULL)
819  debug = true;
820  if (!threeD)
821  nz = 1;
822 
823  /* compute p,q,r from P,Q,R and myid */
824  p = id % px;
825  q = ((id - p) / px) % py;
826  r = (id - p - px * q) / (px * py);
827 
828  if (debug)
829  {
830  sprintf (msgBuf_dh,
831  "this proc's position in subdomain grid: p= %i q= %i r= %i",
832  p, q, r);
833  SET_INFO (msgBuf_dh);
834  }
835 
836  /* compute ilower and iupper from p,q,r and nx,ny,nz */
837  /* zero-based */
838 
839  lowerx = nx * p;
840  upperx = lowerx + nx;
841  lowery = ny * q;
842  uppery = lowery + ny;
843  lowerz = nz * r;
844  upperz = lowerz + nz;
845 
846  if (debug)
847  {
848  sprintf (msgBuf_dh, "local grid parameters: lowerx= %i upperx= %i",
849  lowerx, upperx);
850  SET_INFO (msgBuf_dh);
851  sprintf (msgBuf_dh, "local grid parameters: lowery= %i uppery= %i",
852  lowery, uppery);
853  SET_INFO (msgBuf_dh);
854  sprintf (msgBuf_dh, "local grid parameters: lowerz= %i upperz= %i",
855  lowerz, upperz);
856  SET_INFO (msgBuf_dh);
857  }
858 
859  startRow = mg->first;
860  rp[0] = 0;
861 
862  for (z = lowerz; z < upperz; z++)
863  {
864  for (y = lowery; y < uppery; y++)
865  {
866  for (x = lowerx; x < upperx; x++)
867  {
868 
869  if (debug)
870  {
871  fprintf (logFile, "row= %i x= %i y= %i z= %i\n",
872  localRow + startRow + 1, x, y, z);
873  }
874 
875  /* compute row values and rhs, at the current node */
876  getstencil (mg, x, y, z);
877 
878  /* down plane */
879  if (threeD)
880  {
881  if (z > 0)
882  {
883  naborz1 =
884  rownum (threeD, x, y, z - 1, nx, ny, nz, px, py);
885  cval[idx] = naborz1;
886  aval[idx++] = FRONT (stencil);
887  }
888  }
889 
890  /* south */
891  if (y > 0)
892  {
893  nabory1 = rownum (threeD, x, y - 1, z, nx, ny, nz, px, py);
894  cval[idx] = nabory1;
895  aval[idx++] = SOUTH (stencil);
896  }
897 
898  /* west */
899  if (x > 0)
900  {
901  naborx1 = rownum (threeD, x - 1, y, z, nx, ny, nz, px, py);
902  cval[idx] = naborx1;
903  aval[idx++] = WEST (stencil);
904 /*fprintf(logFile, "--- row: %i; naborx1= %i\n", localRow+startRow+1, 1+naborx1);
905 */
906  }
907 /*
908 else {
909 fprintf(logFile, "--- row: %i; x >= nx*px-1; naborx1 has old value: %i\n", localRow+startRow+1,1+naborx1);
910 }
911 */
912 
913  /* center node */
914  cval[idx] = localRow + startRow;
915  aval[idx++] = CENTER (stencil);
916 
917 
918  /* east */
919  if (x < nx * px - 1)
920  {
921  naborx2 = rownum (threeD, x + 1, y, z, nx, ny, nz, px, py);
922  cval[idx] = naborx2;
923  aval[idx++] = EAST (stencil);
924  }
925 /*
926 else {
927 fprintf(logFile, "--- row: %i; x >= nx*px-1; nobors2 has old value: %i\n", localRow+startRow,1+naborx2);
928 }
929 */
930 
931  /* north */
932  if (y < ny * py - 1)
933  {
934  nabory2 = rownum (threeD, x, y + 1, z, nx, ny, nz, px, py);
935  cval[idx] = nabory2;
936  aval[idx++] = NORTH (stencil);
937  }
938 
939  /* up plane */
940  if (threeD)
941  {
942  if (z < nz * pz - 1)
943  {
944  naborz2 =
945  rownum (threeD, x, y, z + 1, nx, ny, nz, px, py);
946  cval[idx] = naborz2;
947  aval[idx++] = BACK (stencil);
948  }
949  }
950 
951  /* rhs[rhsIdx++] = RHS(stencil); */
952  rhs[localRow] = 0.0;
953 
954  ++localRow;
955  rp[localRow] = idx;
956 
957  /* apply boundary conditions; only for 2D! */
958  if (!threeD && applyBdry)
959  {
960  int globalRow = localRow + startRow - 1;
961  int offset = rp[localRow - 1];
962  int len = rp[localRow] - rp[localRow - 1];
963  double ctr, coeff;
964 
965 /* fprintf(logFile, "globalRow = %i; naborx2 = %i\n", globalRow+1, naborx2+1); */
966 
967  if (x == 0)
968  { /* if x1 */
969  coeff = mg->A (mg->a, x + hhalf, y, z);
970  ctr = mg->A (mg->a, x - hhalf, y, z);
971  setBoundary_private (globalRow, cval + offset,
972  aval + offset, len,
973  &(rhs[localRow - 1]), bcx1, coeff,
974  ctr, naborx2);
975  }
976  else if (x == nx * px - 1)
977  { /* if x2 */
978  coeff = mg->A (mg->a, x - hhalf, y, z);
979  ctr = mg->A (mg->a, x + hhalf, y, z);
980  setBoundary_private (globalRow, cval + offset,
981  aval + offset, len,
982  &(rhs[localRow - 1]), bcx2, coeff,
983  ctr, naborx1);
984  }
985  else if (y == 0)
986  { /* if y1 */
987  coeff = mg->B (mg->b, x, y + hhalf, z);
988  ctr = mg->B (mg->b, x, y - hhalf, z);
989  setBoundary_private (globalRow, cval + offset,
990  aval + offset, len,
991  &(rhs[localRow - 1]), bcy1, coeff,
992  ctr, nabory2);
993  }
994  else if (y == ny * py - 1)
995  { /* if y2 */
996  coeff = mg->B (mg->b, x, y - hhalf, z);
997  ctr = mg->B (mg->b, x, y + hhalf, z);
998  setBoundary_private (globalRow, cval + offset,
999  aval + offset, len,
1000  &(rhs[localRow - 1]), bcy2, coeff,
1001  ctr, nabory1);
1002  }
1003  else if (threeD)
1004  {
1005  if (z == 0)
1006  {
1007  coeff = mg->B (mg->b, x, y, z + hhalf);
1008  ctr = mg->B (mg->b, x, y, z - hhalf);
1009  setBoundary_private (globalRow, cval + offset,
1010  aval + offset, len,
1011  &(rhs[localRow - 1]), bcy1,
1012  coeff, ctr, naborz2);
1013  }
1014  else if (z == nz * nx - 1)
1015  {
1016  coeff = mg->B (mg->b, x, y, z - hhalf);
1017  ctr = mg->B (mg->b, x, y, z + hhalf);
1018  setBoundary_private (globalRow, cval + offset,
1019  aval + offset, len,
1020  &(rhs[localRow - 1]), bcy1,
1021  coeff, ctr, naborz1);
1022  }
1023  }
1024  }
1025  }
1026  }
1027  }
1028 END_FUNC_DH}
1029 
1030 
1031 #undef __FUNC__
1032 #define __FUNC__ "setBoundary_private"
1033 void
1034 setBoundary_private (int node, int *cval, double *aval, int len,
1035  double *rhs, double bc, double coeff, double ctr,
1036  int nabor)
1037 {
1038  START_FUNC_DH int i;
1039 
1040  /* case 1: Dirichlet Boundary condition */
1041  if (bc >= 0)
1042  {
1043  /* set all values to zero, set the diagonal to 1.0, set rhs to "bc" */
1044  *rhs = bc;
1045  for (i = 0; i < len; ++i)
1046  {
1047  if (cval[i] == node)
1048  {
1049  aval[i] = 1.0;
1050  }
1051  else
1052  {
1053  aval[i] = 0;
1054  }
1055  }
1056  }
1057 
1058  /* case 2: neuman */
1059  else
1060  {
1061 /* fprintf(logFile, "node= %i nabor= %i coeff= %g\n", node+1, nabor+1, coeff); */
1062  /* adjust row values */
1063  for (i = 0; i < len; ++i)
1064  {
1065  /* adjust diagonal */
1066  if (cval[i] == node)
1067  {
1068  aval[i] += (ctr - coeff);
1069  /* adust node's right neighbor */
1070  }
1071  else if (cval[i] == nabor)
1072  {
1073  aval[i] = 2.0 * coeff;
1074  }
1075  }
1076  }
1077 END_FUNC_DH}
_mat_dh
Definition: Mat_dh.h:62
_matgenfd
Definition: MatGenFD.h:94
_vec_dh
Definition: Vec_dh.h:52
Generated by   doxygen 1.8.5