MatGenFD.c

00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "MatGenFD.h"
00031 #include "Mat_dh.h"
00032 #include "Vec_dh.h"
00033 #include "Parser_dh.h"
00034 #include "Mem_dh.h"
00035 /* #include "graphColor_dh.h" */
00036 
00037 static bool isThreeD;
00038 
00039   /* handles for values in the 5-point (2D) or 7-point (for 3D) stencil */
00040 #define FRONT(a)  a[5]
00041 #define SOUTH(a)  a[3]
00042 #define WEST(a)   a[1]
00043 #define CENTER(a) a[0]
00044 #define EAST(a)   a[2]
00045 #define NORTH(a)  a[4]
00046 #define BACK(a)   a[6]
00047 #define RHS(a)    a[7]
00048 
00049 static void setBoundary_private (int node, int *cval, double *aval, int len,
00050                  double *rhs, double bc, double coeff,
00051                  double ctr, int nabor);
00052 static void generateStriped (MatGenFD mg, int *rp, int *cval, double *aval,
00053                  Mat_dh A, Vec_dh b);
00054 static void generateBlocked (MatGenFD mg, int *rp, int *cval, double *aval,
00055                  Mat_dh A, Vec_dh b);
00056 static void getstencil (MatGenFD g, int ix, int iy, int iz);
00057 
00058 #if 0
00059 static void fdaddbc (int nx, int ny, int nz, int *rp, int *cval,
00060              int *diag, double *aval, double *rhs, double h,
00061              MatGenFD mg);
00062 #endif
00063 
00064 #undef __FUNC__
00065 #define __FUNC__ "MatGenFDCreate"
00066 void
00067 MatGenFD_Create (MatGenFD * mg)
00068 {
00069   START_FUNC_DH
00070     struct _matgenfd *tmp =
00071     (struct _matgenfd *) MALLOC_DH (sizeof (struct _matgenfd));
00072   CHECK_V_ERROR;
00073   *mg = tmp;
00074 
00075   tmp->debug = Parser_dhHasSwitch (parser_dh, "-debug_matgen");
00076 
00077   tmp->m = 9;
00078   tmp->px = tmp->py = 1;
00079   tmp->pz = 0;
00080   Parser_dhReadInt (parser_dh, "-m", &tmp->m);
00081   Parser_dhReadInt (parser_dh, "-px", &tmp->px);
00082   Parser_dhReadInt (parser_dh, "-py", &tmp->py);
00083   Parser_dhReadInt (parser_dh, "-pz", &tmp->pz);
00084 
00085   if (tmp->px < 1)
00086     tmp->px = 1;
00087   if (tmp->py < 1)
00088     tmp->py = 1;
00089   if (tmp->pz < 0)
00090     tmp->pz = 0;
00091   tmp->threeD = false;
00092   if (tmp->pz)
00093     {
00094       tmp->threeD = true;
00095     }
00096   else
00097     {
00098       tmp->pz = 1;
00099     }
00100   if (Parser_dhHasSwitch (parser_dh, "-threeD"))
00101     tmp->threeD = true;
00102 
00103   tmp->a = tmp->b = tmp->c = 1.0;
00104   tmp->d = tmp->e = tmp->f = 0.0;
00105   tmp->g = tmp->h = 0.0;
00106 
00107   Parser_dhReadDouble (parser_dh, "-dx", &tmp->a);
00108   Parser_dhReadDouble (parser_dh, "-dy", &tmp->b);
00109   Parser_dhReadDouble (parser_dh, "-dz", &tmp->c);
00110   Parser_dhReadDouble (parser_dh, "-cx", &tmp->d);
00111   Parser_dhReadDouble (parser_dh, "-cy", &tmp->e);
00112   Parser_dhReadDouble (parser_dh, "-cz", &tmp->f);
00113 
00114   tmp->a = -1 * fabs (tmp->a);
00115   tmp->b = -1 * fabs (tmp->b);
00116   tmp->c = -1 * fabs (tmp->c);
00117 
00118   tmp->allocateMem = true;
00119 
00120   tmp->A = tmp->B = tmp->C = tmp->D = tmp->E
00121     = tmp->F = tmp->G = tmp->H = konstant;
00122 
00123   tmp->bcX1 = tmp->bcX2 = tmp->bcY1 = tmp->bcY2 = tmp->bcZ1 = tmp->bcZ2 = 0.0;
00124   Parser_dhReadDouble (parser_dh, "-bcx1", &tmp->bcX1);
00125   Parser_dhReadDouble (parser_dh, "-bcx2", &tmp->bcX2);
00126   Parser_dhReadDouble (parser_dh, "-bcy1", &tmp->bcY1);
00127   Parser_dhReadDouble (parser_dh, "-bcy2", &tmp->bcY2);
00128   Parser_dhReadDouble (parser_dh, "-bcz1", &tmp->bcZ1);
00129   Parser_dhReadDouble (parser_dh, "-bcz2", &tmp->bcZ2);
00130 END_FUNC_DH}
00131 
00132 
00133 #undef __FUNC__
00134 #define __FUNC__ "MatGenFD_Destroy"
00135 void
00136 MatGenFD_Destroy (MatGenFD mg)
00137 {
00138   START_FUNC_DH FREE_DH (mg);
00139   CHECK_V_ERROR;
00140 END_FUNC_DH}
00141 
00142 
00143 #undef __FUNC__
00144 #define __FUNC__ "MatGenFD_Run"
00145 void
00146 MatGenFD_Run (MatGenFD mg, int id, int np, Mat_dh * AOut, Vec_dh * rhsOut)
00147 {
00148 /* What this function does:
00149  *   0. creates return objects (A and rhs)
00150  *   1. computes "nice to have" values;
00151  *   2. allocates storage, if required;
00152  *   3. calls generateBlocked() or generateStriped().
00153  *   4. initializes variable in A and rhs.
00154  */
00155 
00156   START_FUNC_DH Mat_dh A;
00157   Vec_dh rhs;
00158   bool threeD = mg->threeD;
00159   int nnz;
00160   int m = mg->m;        /* local unknowns */
00161   bool debug = false, striped;
00162 
00163   if (mg->debug && logFile != NULL)
00164     debug = true;
00165   striped = Parser_dhHasSwitch (parser_dh, "-striped");
00166 
00167   /* 0. create objects */
00168   Mat_dhCreate (AOut);
00169   CHECK_V_ERROR;
00170   Vec_dhCreate (rhsOut);
00171   CHECK_V_ERROR;
00172   A = *AOut;
00173   rhs = *rhsOut;
00174 
00175   /* ensure that processor grid contains the same number of
00176      nodes as there are processors.
00177    */
00178   if (!Parser_dhHasSwitch (parser_dh, "-noChecks"))
00179     {
00180       if (!striped)
00181     {
00182       int npTest = mg->px * mg->py;
00183       if (threeD)
00184         npTest *= mg->pz;
00185       if (npTest != np)
00186         {
00187           sprintf (msgBuf_dh,
00188                "numbers don't match: np_dh = %i, px*py*pz = %i", np,
00189                npTest);
00190           SET_V_ERROR (msgBuf_dh);
00191         }
00192     }
00193     }
00194 
00195   /* 1. compute "nice to have" values */
00196   /* each proc's subgrid dimension */
00197   mg->cc = m;
00198   if (threeD)
00199     {
00200       m = mg->m = m * m * m;
00201     }
00202   else
00203     {
00204       m = mg->m = m * m;
00205     }
00206 
00207   mg->first = id * m;
00208   mg->hh = 1.0 / (mg->px * mg->cc - 1);
00209 
00210   if (debug)
00211     {
00212       sprintf (msgBuf_dh, "cc (local grid dimension) = %i", mg->cc);
00213       SET_INFO (msgBuf_dh);
00214       if (threeD)
00215     {
00216       sprintf (msgBuf_dh, "threeD = true");
00217     }
00218       else
00219     {
00220       sprintf (msgBuf_dh, "threeD = false");
00221     }
00222       SET_INFO (msgBuf_dh);
00223       sprintf (msgBuf_dh, "np= %i  id= %i", np, id);
00224       SET_INFO (msgBuf_dh);
00225     }
00226 
00227   mg->id = id;
00228   mg->np = np;
00229   nnz = threeD ? m * 7 : m * 5;
00230 
00231   /* 2. allocate storage */
00232   if (mg->allocateMem)
00233     {
00234       A->rp = (int *) MALLOC_DH ((m + 1) * sizeof (int));
00235       CHECK_V_ERROR;
00236       A->rp[0] = 0;
00237       A->cval = (int *) MALLOC_DH (nnz * sizeof (int));
00238       CHECK_V_ERROR A->aval = (double *) MALLOC_DH (nnz * sizeof (double));
00239       CHECK_V_ERROR;
00240       /* rhs->vals = (double*)MALLOC_DH(m*sizeof(double)); CHECK_V_ERROR; */
00241     }
00242 
00243   /* 4. initialize variables in A and rhs */
00244   rhs->n = m;
00245   A->m = m;
00246   A->n = m * mg->np;
00247   A->beg_row = mg->first;
00248 
00249   /* 3. generate matrix */
00250   isThreeD = threeD;        /* yuck!  used in box_XX() */
00251   if (Parser_dhHasSwitch (parser_dh, "-striped"))
00252     {
00253       generateStriped (mg, A->rp, A->cval, A->aval, A, rhs);
00254       CHECK_V_ERROR;
00255     }
00256   else
00257     {
00258       generateBlocked (mg, A->rp, A->cval, A->aval, A, rhs);
00259       CHECK_V_ERROR;
00260     }
00261 
00262   /* add in bdry conditions */
00263   /* only implemented for 2D mats! */
00264   if (!threeD)
00265     {
00266 /*  fdaddbc(nx, ny, nz, rp, cval, diag, aval, rhs, h, mg); */
00267     }
00268 
00269 END_FUNC_DH}
00270 
00271 
00272 #undef __FUNC__
00273 #define __FUNC__ "generateStriped"
00274 void
00275 generateStriped (MatGenFD mg, int *rp, int *cval, double *aval, Mat_dh A,
00276          Vec_dh b)
00277 {
00278   START_FUNC_DH int mGlobal;
00279   int m = mg->m;
00280   int beg_row, end_row;
00281   int i, j, k, row;
00282   bool threeD = mg->threeD;
00283   int idx = 0;
00284   double *stencil = mg->stencil;
00285   bool debug = false;
00286   int plane, nodeRemainder;
00287   int naborx1, naborx2, nabory1, nabory2;
00288   double *rhs;
00289 
00290   bool applyBdry = true;
00291   double hhalf;
00292   double bcx1 = mg->bcX1;
00293   double bcx2 = mg->bcX2;
00294   double bcy1 = mg->bcY1;
00295   double bcy2 = mg->bcY2;
00296   /* double bcz1 = mg->bcZ1; */
00297   /* double bcz2 = mg->bcZ2; */
00298   int nx, ny;
00299 
00300   printf_dh ("@@@ using striped partitioning\n");
00301 
00302   if (mg->debug && logFile != NULL)
00303     debug = true;
00304 
00305   /* recompute values (yuck!) */
00306   m = 9;
00307   Parser_dhReadInt (parser_dh, "-m", &m);   /* global grid dimension */
00308   mGlobal = m * m;      /* global unkknowns */
00309   if (threeD)
00310     mGlobal *= m;
00311   i = mGlobal / mg->np;     /* unknowns per processor */
00312   beg_row = i * mg->id;     /* global number of 1st local row */
00313   end_row = beg_row + i;
00314   if (mg->id == mg->np - 1)
00315     end_row = mGlobal;
00316   nx = ny = m;
00317 
00318   mg->hh = 1.0 / (m - 1);
00319   hhalf = 0.5 * mg->hh;
00320 
00321   A->n = m * m;
00322   A->m = end_row - beg_row;
00323   A->beg_row = beg_row;
00324 
00325   Vec_dhInit (b, A->m);
00326   CHECK_V_ERROR;
00327   rhs = b->vals;
00328 
00329   plane = m * m;
00330 
00331   if (debug)
00332     {
00333       fprintf (logFile, "generateStriped: beg_row= %i; end_row= %i; m= %i\n",
00334            beg_row + 1, end_row + 1, m);
00335     }
00336 
00337   for (row = beg_row; row < end_row; ++row)
00338     {
00339       int localRow = row - beg_row;
00340 
00341       /* compute current node's position in grid */
00342       k = (row / plane);
00343       nodeRemainder = row - (k * plane);    /* map row to 1st plane */
00344       j = nodeRemainder / m;
00345       i = nodeRemainder % m;
00346 
00347       if (debug)
00348     {
00349       fprintf (logFile, "row= %i  x= %i  y= %i  z= %i\n", row + 1, i, j,
00350            k);
00351     }
00352 
00353       /* compute column values and rhs entry for the current node */
00354       getstencil (mg, i, j, k);
00355 
00356       /* only homogenous Dirichlet boundary conditions presently supported */
00357 
00358       /* down plane */
00359       if (threeD)
00360     {
00361       if (k > 0)
00362         {
00363           cval[idx] = row - plane;
00364           aval[idx++] = BACK (stencil);
00365         }
00366     }
00367 
00368       /* south */
00369       if (j > 0)
00370     {
00371       nabory1 = cval[idx] = row - m;
00372       aval[idx++] = SOUTH (stencil);
00373     }
00374 
00375       /* west */
00376       if (i > 0)
00377     {
00378       naborx1 = cval[idx] = row - 1;
00379       aval[idx++] = WEST (stencil);
00380     }
00381 
00382       /* center node */
00383       cval[idx] = row;
00384       aval[idx++] = CENTER (stencil);
00385 
00386       /* east */
00387       if (i < m - 1)
00388     {
00389       naborx2 = cval[idx] = row + 1;
00390       aval[idx++] = EAST (stencil);
00391     }
00392 
00393       /* north */
00394       if (j < m - 1)
00395     {
00396       nabory2 = cval[idx] = row + m;
00397       aval[idx++] = NORTH (stencil);
00398     }
00399 
00400       /* up plane */
00401       if (threeD)
00402     {
00403       if (k < m - 1)
00404         {
00405           cval[idx] = row + plane;
00406           aval[idx++] = FRONT (stencil);
00407         }
00408     }
00409       rhs[localRow] = 0.0;
00410       ++localRow;
00411       rp[localRow] = idx;
00412 
00413       /* apply boundary conditions; only for 2D! */
00414       if (!threeD && applyBdry)
00415     {
00416       int offset = rp[localRow - 1];
00417       int len = rp[localRow] - rp[localRow - 1];
00418       double ctr, coeff;
00419 
00420 /* fprintf(logFile, "globalRow = %i; naborx2 = %i\n", row+1, row); */
00421 
00422       if (i == 0)
00423         {           /* if x1 */
00424           coeff = mg->A (mg->a, i + hhalf, j, k);
00425           ctr = mg->A (mg->a, i - hhalf, j, k);
00426           setBoundary_private (row, cval + offset, aval + offset, len,
00427                    &(rhs[localRow - 1]), bcx1, coeff, ctr,
00428                    naborx2);
00429         }
00430       else if (i == nx - 1)
00431         {           /* if x2 */
00432           coeff = mg->A (mg->a, i - hhalf, j, k);
00433           ctr = mg->A (mg->a, i + hhalf, j, k);
00434           setBoundary_private (row, cval + offset, aval + offset, len,
00435                    &(rhs[localRow - 1]), bcx2, coeff, ctr,
00436                    naborx1);
00437         }
00438       else if (j == 0)
00439         {           /* if y1 */
00440           coeff = mg->B (mg->b, i, j + hhalf, k);
00441           ctr = mg->B (mg->b, i, j - hhalf, k);
00442           setBoundary_private (row, cval + offset, aval + offset, len,
00443                    &(rhs[localRow - 1]), bcy1, coeff, ctr,
00444                    nabory2);
00445         }
00446       else if (j == ny - 1)
00447         {           /* if y2 */
00448           coeff = mg->B (mg->b, i, j - hhalf, k);
00449           ctr = mg->B (mg->b, i, j + hhalf, k);
00450           setBoundary_private (row, cval + offset, aval + offset, len,
00451                    &(rhs[localRow - 1]), bcy2, coeff, ctr,
00452                    nabory1);
00453         }
00454     }
00455     }
00456 END_FUNC_DH}
00457 
00458 
00459 /* zero-based 
00460    (from Edmond Chow)
00461 */
00462 /* 
00463    x,y,z       -  coordinates of row, wrt naturally ordered grid
00464    nz, ny, nz  -  local grid dimensions, wrt 0
00465    P, Q        -  subdomain grid dimensions in x and y directions
00466 */
00467 int
00468 rownum (const bool threeD, const int x, const int y, const int z,
00469     const int nx, const int ny, const int nz, int P, int Q)
00470 {
00471   int p, q, r;
00472   int lowerx, lowery, lowerz;
00473   int id, startrow;
00474 
00475 
00476   /* compute x,y,z coordinates of subdomain to which
00477      this row belongs.
00478    */
00479   p = x / nx;
00480   q = y / ny;
00481   r = z / nz;
00482 
00483 /*
00484 if (myid_dh == 0) printf("nx= %i  ny= %i  nz= %i\n", nx, ny, nz);
00485 if (myid_dh == 0) printf("x= %i y= %i z= %i  threeD= %i  p= %i q= %i r= %i\n",
00486               x,y,z,threeD, p,q,r);
00487 */
00488 
00489   /* compute the subdomain (processor) of the subdomain to which
00490      this row belongs.
00491    */
00492   if (threeD)
00493     {
00494       id = r * P * Q + q * P + p;
00495     }
00496   else
00497     {
00498       id = q * P + p;
00499     }
00500 
00501 /*  if (myid_dh == 0) printf(" id= %i\n", id);
00502 */
00503 
00504   /* smallest row in the subdomain */
00505   startrow = id * (nx * ny * nz);
00506 
00507   /* x,y, and z coordinates of local grid of unknowns */
00508   lowerx = nx * p;
00509   lowery = ny * q;
00510   lowerz = nz * r;
00511 
00512   if (threeD)
00513     {
00514       return startrow + nx * ny * (z - lowerz) + nx * (y - lowery) + (x -
00515                                       lowerx);
00516     }
00517   else
00518     {
00519       return startrow + nx * (y - lowery) + (x - lowerx);
00520     }
00521 }
00522 
00523 
00524 
00525 void
00526 getstencil (MatGenFD g, int ix, int iy, int iz)
00527 {
00528   int k;
00529   double h = g->hh;
00530   double hhalf = h * 0.5;
00531   double x = h * ix;
00532   double y = h * iy;
00533   double z = h * iz;
00534   double cntr = 0.0;
00535   double *stencil = g->stencil;
00536   double coeff;
00537   bool threeD = g->threeD;
00538 
00539   for (k = 0; k < 8; ++k)
00540     stencil[k] = 0.0;
00541 
00542   /* differentiation wrt x */
00543   coeff = g->A (g->a, x + hhalf, y, z);
00544   EAST (stencil) += coeff;
00545   cntr += coeff;
00546 
00547   coeff = g->A (g->a, x - hhalf, y, z);
00548   WEST (stencil) += coeff;
00549   cntr += coeff;
00550 
00551   coeff = g->D (g->d, x, y, z) * hhalf;
00552   EAST (stencil) += coeff;
00553   WEST (stencil) -= coeff;
00554 
00555   /* differentiation wrt y */
00556   coeff = g->B (g->b, x, y + hhalf, z);
00557   NORTH (stencil) += coeff;
00558   cntr += coeff;
00559 
00560   coeff = g->B (g->b, x, y - hhalf, z);
00561   SOUTH (stencil) += coeff;
00562   cntr += coeff;
00563 
00564   coeff = g->E (g->e, x, y, z) * hhalf;
00565   NORTH (stencil) += coeff;
00566   SOUTH (stencil) -= coeff;
00567 
00568   /* differentiation wrt z */
00569   if (threeD)
00570     {
00571       coeff = g->C (g->c, x, y, z + hhalf);
00572       BACK (stencil) += coeff;
00573       cntr += coeff;
00574 
00575       coeff = g->C (g->c, x, y, z - hhalf);
00576       FRONT (stencil) += coeff;
00577       cntr += coeff;
00578 
00579       coeff = g->F (g->f, x, y, z) * hhalf;
00580       BACK (stencil) += coeff;
00581       FRONT (stencil) -= coeff;
00582     }
00583 
00584   /* contribution from function G: */
00585   coeff = g->G (g->g, x, y, z);
00586   CENTER (stencil) = h * h * coeff - cntr;
00587 
00588   RHS (stencil) = h * h * g->H (g->h, x, y, z);
00589 }
00590 
00591 
00592 double
00593 konstant (double coeff, double x, double y, double z)
00594 {
00595   return coeff;
00596 }
00597 
00598 double
00599 e2_xy (double coeff, double x, double y, double z)
00600 {
00601   return exp (coeff * x * y);
00602 }
00603 
00604 double boxThreeD (double coeff, double x, double y, double z);
00605 
00606 /* returns diffusivity constant -bd1 if the point
00607    (x,y,z) is inside the box whose upper left and
00608    lower right points are (-bx1,-by1), (-bx2,-by2);
00609    else, returns diffusivity constant -bd2
00610 */
00611 double
00612 box_1 (double coeff, double x, double y, double z)
00613 {
00614   static bool setup = false;
00615   double retval = coeff;
00616 
00617   /* dffusivity constants */
00618   static double dd1 = BOX1_DD;
00619   static double dd2 = BOX2_DD;
00620   static double dd3 = BOX3_DD;
00621 
00622   /* boxes */
00623   static double ax1 = BOX1_X1, ay1 = BOX1_Y1;
00624   static double ax2 = BOX1_X2, ay2 = BOX1_Y2;
00625   static double bx1 = BOX2_X1, by1 = BOX2_Y1;
00626   static double bx2 = BOX2_X2, by2 = BOX2_Y2;
00627   static double cx1 = BOX3_X1, cy1 = BOX3_Y1;
00628   static double cx2 = BOX3_X2, cy2 = BOX3_Y2;
00629 
00630   if (isThreeD)
00631     {
00632       return (boxThreeD (coeff, x, y, z));
00633     }
00634 
00635 
00636   /* 1st time through, parse for dffusivity constants */
00637   if (!setup)
00638     {
00639       dd1 = 0.1;
00640       dd2 = 0.1;
00641       dd3 = 10;
00642       Parser_dhReadDouble (parser_dh, "-dd1", &dd1);
00643       Parser_dhReadDouble (parser_dh, "-dd2", &dd2);
00644       Parser_dhReadDouble (parser_dh, "-dd3", &dd3);
00645       Parser_dhReadDouble (parser_dh, "-box1x1", &cx1);
00646       Parser_dhReadDouble (parser_dh, "-box1x2", &cx2);
00647       setup = true;
00648     }
00649 
00650   /* determine if point is inside box a */
00651   if (x > ax1 && x < ax2 && y > ay1 && y < ay2)
00652     {
00653       retval = dd1 * coeff;
00654     }
00655 
00656   /* determine if point is inside box b */
00657   if (x > bx1 && x < bx2 && y > by1 && y < by2)
00658     {
00659       retval = dd2 * coeff;
00660     }
00661 
00662   /* determine if point is inside box c */
00663   if (x > cx1 && x < cx2 && y > cy1 && y < cy2)
00664     {
00665       retval = dd3 * coeff;
00666     }
00667 
00668   return retval;
00669 }
00670 
00671 double
00672 boxThreeD (double coeff, double x, double y, double z)
00673 {
00674   static bool setup = false;
00675   double retval = coeff;
00676 
00677   /* dffusivity constants */
00678   static double dd1 = 100;
00679 
00680   /* boxes */
00681   static double x1 = .2, x2 = .8;
00682   static double y1 = .3, y2 = .7;
00683   static double z1 = .4, z2 = .6;
00684 
00685   /* 1st time through, parse for diffusivity constants */
00686   if (!setup)
00687     {
00688       Parser_dhReadDouble (parser_dh, "-dd1", &dd1);
00689       setup = true;
00690     }
00691 
00692   /* determine if point is inside the box */
00693   if (x > x1 && x < x2 && y > y1 && y < y2 && z > z1 && z < z2)
00694     {
00695       retval = dd1 * coeff;
00696     }
00697 
00698   return retval;
00699 }
00700 
00701 #if 0
00702 double
00703 box_1 (double coeff, double x, double y, double z)
00704 {
00705   static double x1, x2, y1, y2;
00706   static double d1, d2;
00707   bool setup = false;
00708   double retval;
00709 
00710   /* 1st time through, parse for constants and
00711      bounding box definition
00712    */
00713   if (!setup)
00714     {
00715       x1 = .25;
00716       x2 = .75;
00717       y1 = .25;
00718       y2 = .75;
00719       d1 = 1;
00720       d2 = 2;
00721       Parser_dhReadDouble (parser_dh, "-bx1", &x1);
00722       Parser_dhReadDouble (parser_dh, "-bx2", &x2);
00723       Parser_dhReadDouble (parser_dh, "-by1", &y1);
00724       Parser_dhReadDouble (parser_dh, "-by2", &y2);
00725       Parser_dhReadDouble (parser_dh, "-bd1", &d1);
00726       Parser_dhReadDouble (parser_dh, "-bd2", &d2);
00727       setup = true;
00728     }
00729 
00730   retval = d2;
00731 
00732   /* determine if point is inside box */
00733   if (x > x1 && x < x2 && y > y1 && y < y2)
00734     {
00735       retval = d1;
00736     }
00737 
00738   return -1 * retval;
00739 }
00740 #endif
00741 
00742 /* divide square into 4 quadrants; return one of
00743    2 constants depending on the quadrant (checkerboard)
00744 */
00745 double
00746 box_2 (double coeff, double x, double y, double z)
00747 {
00748   bool setup = false;
00749   static double d1, d2;
00750   double retval;
00751 
00752   if (!setup)
00753     {
00754       d1 = 1;
00755       d2 = 2;
00756       Parser_dhReadDouble (parser_dh, "-bd1", &d1);
00757       Parser_dhReadDouble (parser_dh, "-bd2", &d2);
00758     }
00759 
00760   retval = d2;
00761 
00762   if (x < .5 && y < .5)
00763     retval = d1;
00764   if (x > .5 && y > .5)
00765     retval = d1;
00766 
00767   return -1 * retval;
00768 }
00769 
00770 
00771 #undef __FUNC__
00772 #define __FUNC__ "generateBlocked"
00773 void
00774 generateBlocked (MatGenFD mg, int *rp, int *cval, double *aval, Mat_dh A,
00775          Vec_dh b)
00776 {
00777   START_FUNC_DH bool applyBdry = true;
00778   double *stencil = mg->stencil;
00779   int id = mg->id;
00780   bool threeD = mg->threeD;
00781   int px = mg->px, py = mg->py, pz = mg->pz;    /* processor grid dimensions */
00782   int p, q, r;          /* this proc's position in processor grid */
00783   int cc = mg->cc;      /* local grid dimension (grid of unknowns) */
00784   int nx = cc, ny = cc, nz = cc;
00785   int lowerx, upperx, lowery, uppery, lowerz, upperz;
00786   int startRow;
00787   int x, y, z;
00788   bool debug = false;
00789   int idx = 0, localRow = 0;    /* nabor; */
00790   int naborx1, naborx2, nabory1, nabory2, naborz1, naborz2;
00791   double *rhs;
00792 
00793   double hhalf = 0.5 * mg->hh;
00794   double bcx1 = mg->bcX1;
00795   double bcx2 = mg->bcX2;
00796   double bcy1 = mg->bcY1;
00797   double bcy2 = mg->bcY2;
00798   /* double bcz1 = mg->bcZ1; */
00799   /* double bcz2 = mg->bcZ2; */
00800 
00801   Vec_dhInit (b, A->m);
00802   CHECK_V_ERROR;
00803   rhs = b->vals;
00804 
00805   if (mg->debug && logFile != NULL)
00806     debug = true;
00807   if (!threeD)
00808     nz = 1;
00809 
00810   /* compute p,q,r from P,Q,R and myid */
00811   p = id % px;
00812   q = ((id - p) / px) % py;
00813   r = (id - p - px * q) / (px * py);
00814 
00815   if (debug)
00816     {
00817       sprintf (msgBuf_dh,
00818            "this proc's position in subdomain grid: p= %i  q= %i  r= %i",
00819            p, q, r);
00820       SET_INFO (msgBuf_dh);
00821     }
00822 
00823   /* compute ilower and iupper from p,q,r and nx,ny,nz */
00824   /* zero-based */
00825 
00826   lowerx = nx * p;
00827   upperx = lowerx + nx;
00828   lowery = ny * q;
00829   uppery = lowery + ny;
00830   lowerz = nz * r;
00831   upperz = lowerz + nz;
00832 
00833   if (debug)
00834     {
00835       sprintf (msgBuf_dh, "local grid parameters: lowerx= %i  upperx= %i",
00836            lowerx, upperx);
00837       SET_INFO (msgBuf_dh);
00838       sprintf (msgBuf_dh, "local grid parameters: lowery= %i  uppery= %i",
00839            lowery, uppery);
00840       SET_INFO (msgBuf_dh);
00841       sprintf (msgBuf_dh, "local grid parameters: lowerz= %i  upperz= %i",
00842            lowerz, upperz);
00843       SET_INFO (msgBuf_dh);
00844     }
00845 
00846   startRow = mg->first;
00847   rp[0] = 0;
00848 
00849   for (z = lowerz; z < upperz; z++)
00850     {
00851       for (y = lowery; y < uppery; y++)
00852     {
00853       for (x = lowerx; x < upperx; x++)
00854         {
00855 
00856           if (debug)
00857         {
00858           fprintf (logFile, "row= %i  x= %i  y= %i  z= %i\n",
00859                localRow + startRow + 1, x, y, z);
00860         }
00861 
00862           /* compute row values and rhs, at the current node */
00863           getstencil (mg, x, y, z);
00864 
00865           /* down plane */
00866           if (threeD)
00867         {
00868           if (z > 0)
00869             {
00870               naborz1 =
00871             rownum (threeD, x, y, z - 1, nx, ny, nz, px, py);
00872               cval[idx] = naborz1;
00873               aval[idx++] = FRONT (stencil);
00874             }
00875         }
00876 
00877           /* south */
00878           if (y > 0)
00879         {
00880           nabory1 = rownum (threeD, x, y - 1, z, nx, ny, nz, px, py);
00881           cval[idx] = nabory1;
00882           aval[idx++] = SOUTH (stencil);
00883         }
00884 
00885           /* west */
00886           if (x > 0)
00887         {
00888           naborx1 = rownum (threeD, x - 1, y, z, nx, ny, nz, px, py);
00889           cval[idx] = naborx1;
00890           aval[idx++] = WEST (stencil);
00891 /*fprintf(logFile, "--- row: %i;  naborx1= %i\n", localRow+startRow+1, 1+naborx1);
00892 */
00893         }
00894 /*
00895 else {
00896 fprintf(logFile, "--- row: %i;  x >= nx*px-1; naborx1 has old value: %i\n", localRow+startRow+1,1+naborx1);
00897 }
00898 */
00899 
00900           /* center node */
00901           cval[idx] = localRow + startRow;
00902           aval[idx++] = CENTER (stencil);
00903 
00904 
00905           /* east */
00906           if (x < nx * px - 1)
00907         {
00908           naborx2 = rownum (threeD, x + 1, y, z, nx, ny, nz, px, py);
00909           cval[idx] = naborx2;
00910           aval[idx++] = EAST (stencil);
00911         }
00912 /*
00913 else {
00914 fprintf(logFile, "--- row: %i;  x >= nx*px-1; nobors2 has old value: %i\n", localRow+startRow,1+naborx2);
00915 }
00916 */
00917 
00918           /* north */
00919           if (y < ny * py - 1)
00920         {
00921           nabory2 = rownum (threeD, x, y + 1, z, nx, ny, nz, px, py);
00922           cval[idx] = nabory2;
00923           aval[idx++] = NORTH (stencil);
00924         }
00925 
00926           /* up plane */
00927           if (threeD)
00928         {
00929           if (z < nz * pz - 1)
00930             {
00931               naborz2 =
00932             rownum (threeD, x, y, z + 1, nx, ny, nz, px, py);
00933               cval[idx] = naborz2;
00934               aval[idx++] = BACK (stencil);
00935             }
00936         }
00937 
00938           /* rhs[rhsIdx++] = RHS(stencil); */
00939           rhs[localRow] = 0.0;
00940 
00941           ++localRow;
00942           rp[localRow] = idx;
00943 
00944           /* apply boundary conditions; only for 2D! */
00945           if (!threeD && applyBdry)
00946         {
00947           int globalRow = localRow + startRow - 1;
00948           int offset = rp[localRow - 1];
00949           int len = rp[localRow] - rp[localRow - 1];
00950           double ctr, coeff;
00951 
00952 /* fprintf(logFile, "globalRow = %i; naborx2 = %i\n", globalRow+1, naborx2+1); */
00953 
00954           if (x == 0)
00955             {       /* if x1 */
00956               coeff = mg->A (mg->a, x + hhalf, y, z);
00957               ctr = mg->A (mg->a, x - hhalf, y, z);
00958               setBoundary_private (globalRow, cval + offset,
00959                        aval + offset, len,
00960                        &(rhs[localRow - 1]), bcx1, coeff,
00961                        ctr, naborx2);
00962             }
00963           else if (x == nx * px - 1)
00964             {       /* if x2 */
00965               coeff = mg->A (mg->a, x - hhalf, y, z);
00966               ctr = mg->A (mg->a, x + hhalf, y, z);
00967               setBoundary_private (globalRow, cval + offset,
00968                        aval + offset, len,
00969                        &(rhs[localRow - 1]), bcx2, coeff,
00970                        ctr, naborx1);
00971             }
00972           else if (y == 0)
00973             {       /* if y1 */
00974               coeff = mg->B (mg->b, x, y + hhalf, z);
00975               ctr = mg->B (mg->b, x, y - hhalf, z);
00976               setBoundary_private (globalRow, cval + offset,
00977                        aval + offset, len,
00978                        &(rhs[localRow - 1]), bcy1, coeff,
00979                        ctr, nabory2);
00980             }
00981           else if (y == ny * py - 1)
00982             {       /* if y2 */
00983               coeff = mg->B (mg->b, x, y - hhalf, z);
00984               ctr = mg->B (mg->b, x, y + hhalf, z);
00985               setBoundary_private (globalRow, cval + offset,
00986                        aval + offset, len,
00987                        &(rhs[localRow - 1]), bcy2, coeff,
00988                        ctr, nabory1);
00989             }
00990           else if (threeD)
00991             {
00992               if (z == 0)
00993             {
00994               coeff = mg->B (mg->b, x, y, z + hhalf);
00995               ctr = mg->B (mg->b, x, y, z - hhalf);
00996               setBoundary_private (globalRow, cval + offset,
00997                            aval + offset, len,
00998                            &(rhs[localRow - 1]), bcy1,
00999                            coeff, ctr, naborz2);
01000             }
01001               else if (z == nz * nx - 1)
01002             {
01003               coeff = mg->B (mg->b, x, y, z - hhalf);
01004               ctr = mg->B (mg->b, x, y, z + hhalf);
01005               setBoundary_private (globalRow, cval + offset,
01006                            aval + offset, len,
01007                            &(rhs[localRow - 1]), bcy1,
01008                            coeff, ctr, naborz1);
01009             }
01010             }
01011         }
01012         }
01013     }
01014     }
01015 END_FUNC_DH}
01016 
01017 
01018 #undef __FUNC__
01019 #define __FUNC__ "setBoundary_private"
01020 void
01021 setBoundary_private (int node, int *cval, double *aval, int len,
01022              double *rhs, double bc, double coeff, double ctr,
01023              int nabor)
01024 {
01025   START_FUNC_DH int i;
01026 
01027   /* case 1: Dirichlet Boundary condition  */
01028   if (bc >= 0)
01029     {
01030       /* set all values to zero, set the diagonal to 1.0, set rhs to "bc" */
01031       *rhs = bc;
01032       for (i = 0; i < len; ++i)
01033     {
01034       if (cval[i] == node)
01035         {
01036           aval[i] = 1.0;
01037         }
01038       else
01039         {
01040           aval[i] = 0;
01041         }
01042     }
01043     }
01044 
01045   /* case 2: neuman */
01046   else
01047     {
01048 /* fprintf(logFile, "node= %i  nabor= %i  coeff= %g\n", node+1, nabor+1, coeff); */
01049       /* adjust row values */
01050       for (i = 0; i < len; ++i)
01051     {
01052       /* adjust diagonal */
01053       if (cval[i] == node)
01054         {
01055           aval[i] += (ctr - coeff);
01056           /* adust node's right neighbor */
01057         }
01058       else if (cval[i] == nabor)
01059         {
01060           aval[i] = 2.0 * coeff;
01061         }
01062     }
01063     }
01064 END_FUNC_DH}

Generated on Tue Jul 13 09:27:14 2010 for IFPACK by  doxygen 1.4.7