IFPACK Development
Loading...
Searching...
No Matches
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
50static 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
62static void setBoundary_private (int node, int *cval, double *aval, int len,
63 double *rhs, double bc, double coeff,
64 double ctr, int nabor);
65static void generateStriped (MatGenFD mg, int *rp, int *cval, double *aval,
66 Mat_dh A, Vec_dh b);
67static void generateBlocked (MatGenFD mg, int *rp, int *cval, double *aval,
68 Mat_dh A, Vec_dh b);
69static void getstencil (MatGenFD g, int ix, int iy, int iz);
70
71#if 0
72static 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"
79void
80MatGenFD_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);
143END_FUNC_DH}
144
145
146#undef __FUNC__
147#define __FUNC__ "MatGenFD_Destroy"
148void
149MatGenFD_Destroy (MatGenFD mg)
150{
151 START_FUNC_DH FREE_DH (mg);
152 CHECK_V_ERROR;
153END_FUNC_DH}
154
155
156#undef __FUNC__
157#define __FUNC__ "MatGenFD_Run"
158void
159MatGenFD_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
282END_FUNC_DH}
283
284
285#undef __FUNC__
286#define __FUNC__ "generateStriped"
287void
288generateStriped (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 }
469END_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*/
480int
481rownum (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/*
497if (myid_dh == 0) printf("nx= %i ny= %i nz= %i\n", nx, ny, nz);
498if (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
538void
539getstencil (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
605double
606konstant (double coeff, double x, double y, double z)
607{
608 return coeff;
609}
610
611double
612e2_xy (double coeff, double x, double y, double z)
613{
614 return exp (coeff * x * y);
615}
616
617double 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*/
624double
625box_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
684double
685boxThreeD (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
715double
716box_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*/
758double
759box_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"
786void
787generateBlocked (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/*
908else {
909fprintf(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/*
926else {
927fprintf(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 }
1028END_FUNC_DH}
1029
1030
1031#undef __FUNC__
1032#define __FUNC__ "setBoundary_private"
1033void
1034setBoundary_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 }
1077END_FUNC_DH}
Definition: Mat_dh.h:63
Definition: Vec_dh.h:53