Ifpack Package Browser (Single Doxygen Collection) Development
MatGenFD.c
Go to the documentation of this file.
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}
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines