Actual source code: ex2.c

      1: static char help[] = "This example is based on ex1.c, but generates a random network of chosen sizes and parameters. \n\
      2:   Usage: -n determines number of nodes. The nonnegative seed can be specified with the flag -seed, otherwise the program generates a random seed.\n\n";

      4: /* T
      5:   Concepts: DMNetwork
      6:   Concepts: KSP
      7: */

      9: #include <petscdmnetwork.h>
     10: #include <petscksp.h>
     11: #include <time.h>

     13: typedef struct {
     14:   PetscInt    id; /* node id */
     15:   PetscScalar inj; /* current injection (A) */
     16:   PetscBool   gr; /* grounded ? */
     17: } Node;

     19: typedef struct {
     20:   PetscInt    id;  /* branch id */
     21:   PetscScalar r;   /* resistance (ohms) */
     22:   PetscScalar bat; /* battery (V) */
     23: } Branch;

     25: typedef struct Edge {
     26:   PetscInt      n;
     27:   PetscInt      i;
     28:   PetscInt      j;
     29:   struct Edge   *next;
     30: } Edge;

     32: PetscReal findDistance(PetscReal x1, PetscReal x2, PetscReal y1, PetscReal y2)
     33: {
     34:   return PetscSqrtReal(PetscPowReal(x2-x1,2.0) + PetscPowReal(y2-y1,2.0));
     35: }

     37: /*
     38:   The algorithm for network formation is based on the paper:
     39:   Routing of Multipoint Connections, Bernard M. Waxman. 1988
     40: */

     42: PetscErrorCode random_network(PetscInt nvertex,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist,PetscInt seed)
     43: {
     45:   PetscInt       i, j, nedges = 0;
     46:   PetscInt       *edgelist;
     47:   PetscInt       nbat, ncurr, fr, to;
     48:   PetscReal      *x, *y, value, xmax = 10.0; /* generate points in square */
     49:   PetscReal      maxdist = 0.0, dist, alpha, beta, prob;
     50:   PetscRandom    rnd;
     51:   Branch         *branch;
     52:   Node           *node;
     53:   Edge           *head = NULL, *nnew= NULL, *aux= NULL;

     56:   PetscRandomCreate(PETSC_COMM_SELF,&rnd);
     57:   PetscRandomSetFromOptions(rnd);

     59:   PetscRandomSetSeed(rnd, seed);
     60:   PetscRandomSeed(rnd);

     62:   /* These parameters might be modified for experimentation */
     63:   nbat  = (PetscInt)(0.1*nvertex);
     64:   ncurr = (PetscInt)(0.1*nvertex);
     65:   alpha = 0.6;
     66:   beta  = 0.2;

     68:   PetscMalloc2(nvertex,&x,nvertex,&y);

     70:   PetscRandomSetInterval(rnd,0.0,xmax);
     71:   for (i=0; i<nvertex; i++) {
     72:     PetscRandomGetValueReal(rnd,&x[i]);
     73:     PetscRandomGetValueReal(rnd,&y[i]);
     74:   }

     76:   /* find maximum distance */
     77:   for (i=0; i<nvertex; i++) {
     78:     for (j=0; j<nvertex; j++) {
     79:       dist = findDistance(x[i],x[j],y[i],y[j]);
     80:       if (dist >= maxdist) maxdist = dist;
     81:     }
     82:   }

     84:   PetscRandomSetInterval(rnd,0.0,1.0);
     85:   for (i=0; i<nvertex; i++) {
     86:     for (j=0; j<nvertex; j++) {
     87:       if (j != i) {
     88:         dist = findDistance(x[i],x[j],y[i],y[j]);
     89:         prob = beta*PetscExpScalar(-dist/(maxdist*alpha));
     90:         PetscRandomGetValueReal(rnd,&value);
     91:         if (value <= prob) {
     92:           PetscMalloc1(1,&nnew);
     93:           if (head == NULL) {
     94:             head       = nnew;
     95:             head->next = NULL;
     96:             head->n    = nedges;
     97:             head->i    = i;
     98:             head->j    = j;
     99:           } else {
    100:             aux = head;
    101:             head = nnew;
    102:             head->n    = nedges;
    103:             head->next = aux;
    104:             head->i    = i;
    105:             head->j    = j;
    106:           }
    107:           nedges += 1;
    108:         }
    109:       }
    110:     }
    111:   }

    113:   PetscMalloc1(2*nedges,&edgelist);

    115:   for (aux = head; aux; aux = aux->next) {
    116:     edgelist[(aux->n)*2]     = aux->i;
    117:     edgelist[(aux->n)*2 + 1] = aux->j;
    118:   }

    120:   aux = head;
    121:   while (aux != NULL) {
    122:     nnew = aux;
    123:     aux = aux->next;
    124:     PetscFree(nnew);
    125:   }

    127:   PetscCalloc2(nvertex,&node,nedges,&branch);

    129:   for (i = 0; i < nvertex; i++) {
    130:     node[i].id  = i;
    131:     node[i].inj = 0;
    132:     node[i].gr = PETSC_FALSE;
    133:   }

    135:   for (i = 0; i < nedges; i++) {
    136:     branch[i].id  = i;
    137:     branch[i].r   = 1.0;
    138:     branch[i].bat = 0;
    139:   }

    141:   /* Chose random node as ground voltage */
    142:   PetscRandomSetInterval(rnd,0.0,nvertex);
    143:   PetscRandomGetValueReal(rnd,&value);
    144:   node[(int)value].gr = PETSC_TRUE;

    146:   /* Create random current and battery injectionsa */
    147:   for (i=0; i<ncurr; i++) {
    148:     PetscRandomSetInterval(rnd,0.0,nvertex);
    149:     PetscRandomGetValueReal(rnd,&value);
    150:     fr   = edgelist[(int)value*2];
    151:     to   = edgelist[(int)value*2 + 1];
    152:     node[fr].inj += 1.0;
    153:     node[to].inj -= 1.0;
    154:   }

    156:   for (i=0; i<nbat; i++) {
    157:     PetscRandomSetInterval(rnd,0.0,nedges);
    158:     PetscRandomGetValueReal(rnd,&value);
    159:     branch[(int)value].bat += 1.0;
    160:   }

    162:   PetscFree2(x,y);
    163:   PetscRandomDestroy(&rnd);

    165:   /* assign pointers */
    166:   *pnbranch  = nedges;
    167:   *pedgelist = edgelist;
    168:   *pbranch   = branch;
    169:   *pnode     = node;
    170:   PetscFunctionReturn(ierr);
    171: }

    173: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
    174: {
    175:   PetscErrorCode    ierr;
    176:   Vec               localb;
    177:   Branch            *branch;
    178:   Node              *node;
    179:   PetscInt          e,v,vStart,vEnd,eStart, eEnd;
    180:   PetscInt          lofst,lofst_to,lofst_fr,row[2],col[6];
    181:   PetscBool         ghost;
    182:   const PetscInt    *cone;
    183:   PetscScalar       *barr,val[6];

    186:   DMGetLocalVector(networkdm,&localb);
    187:   VecSet(b,0.0);
    188:   VecSet(localb,0.0);
    189:   MatZeroEntries(A);

    191:   VecGetArray(localb,&barr);

    193:   /*
    194:     We can define the current as a "edge characteristic" and the voltage
    195:     and the voltage as a "vertex characteristic". With that, we can iterate
    196:     the list of edges and vertices, query the associated voltages and currents
    197:     and use them to write the Kirchoff equations.
    198:   */

    200:   /* Branch equations: i/r + uj - ui = battery */
    201:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    202:   for (e = 0; e < eEnd; e++) {
    203:     DMNetworkGetComponent(networkdm,e,0,NULL,(void**)&branch,NULL);
    204:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);

    206:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
    207:     DMNetworkGetLocalVecOffset(networkdm,cone[0],ALL_COMPONENTS,&lofst_fr);
    208:     DMNetworkGetLocalVecOffset(networkdm,cone[1],ALL_COMPONENTS,&lofst_to);

    210:     barr[lofst] = branch->bat;

    212:     row[0] = lofst;
    213:     col[0] = lofst;     val[0] =  1;
    214:     col[1] = lofst_to;  val[1] =  1;
    215:     col[2] = lofst_fr;  val[2] = -1;
    216:     MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);

    218:     /* from node */
    219:     DMNetworkGetComponent(networkdm,cone[0],0,NULL,(void**)&node,NULL);

    221:     if (!node->gr) {
    222:       row[0] = lofst_fr;
    223:       col[0] = lofst;   val[0] =  1;
    224:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
    225:     }

    227:     /* to node */
    228:     DMNetworkGetComponent(networkdm,cone[1],0,NULL,(void**)&node,NULL);

    230:     if (!node->gr) {
    231:       row[0] = lofst_to;
    232:       col[0] = lofst;   val[0] =  -1;
    233:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
    234:     }
    235:   }

    237:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    238:   for (v = vStart; v < vEnd; v++) {
    239:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
    240:     if (!ghost) {
    241:       DMNetworkGetComponent(networkdm,v,0,NULL,(void**)&node,NULL);
    242:       DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);

    244:       if (node->gr) {
    245:         row[0] = lofst;
    246:         col[0] = lofst;   val[0] =  1;
    247:         MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
    248:       } else {
    249:         barr[lofst] -= node->inj;
    250:       }
    251:     }
    252:   }

    254:   VecRestoreArray(localb,&barr);

    256:   DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
    257:   DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
    258:   DMRestoreLocalVector(networkdm,&localb);

    260:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    261:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    262:   return(0);
    263: }

    265: int main(int argc,char ** argv)
    266: {
    267:   PetscErrorCode    ierr;
    268:   PetscInt          i, nbranch = 0, eStart, eEnd, vStart, vEnd;
    269:   PetscInt          seed = 0, nnode = 0;
    270:   PetscMPIInt       size, rank;
    271:   DM                networkdm;
    272:   Vec               x, b;
    273:   Mat               A;
    274:   KSP               ksp;
    275:   PetscInt          *edgelist = NULL;
    276:   PetscInt          componentkey[2];
    277:   Node              *node;
    278:   Branch            *branch;
    279: #if defined(PETSC_USE_LOG)
    280:   PetscLogStage stage[3];
    281: #endif

    283:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
    284:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    285:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

    287:   PetscOptionsGetInt(NULL,NULL,"-seed",&seed,NULL);

    289:   PetscLogStageRegister("Network Creation", &stage[0]);
    290:   PetscLogStageRegister("DMNetwork data structures", &stage[1]);
    291:   PetscLogStageRegister("KSP", &stage[2]);

    293:   PetscLogStagePush(stage[0]);
    294:   /* "read" data only for processor 0 */
    295:   if (!rank) {
    296:     nnode = 100;
    297:     PetscOptionsGetInt(NULL,NULL,"-n",&nnode,NULL);
    298:     random_network(nnode, &nbranch, &node, &branch, &edgelist, seed);
    299:   }
    300:   PetscLogStagePop();

    302:   PetscLogStagePush(stage[1]);
    303:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
    304:   DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
    305:   DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);

    307:   /* Set number of nodes/edges and edge connectivity */
    308:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
    309:   DMNetworkAddSubnetwork(networkdm,"",nnode,nbranch,edgelist,NULL);

    311:   /* Set up the network layout */
    312:   DMNetworkLayoutSetUp(networkdm);

    314:   /* Add network components (physical parameters of nodes and branches) and num of variables */
    315:   if (!rank) {
    316:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    317:     for (i = eStart; i < eEnd; i++) {
    318:       DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
    319:     }

    321:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    322:     for (i = vStart; i < vEnd; i++) {
    323:       DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
    324:     }
    325:   }

    327:   /* Network partitioning and distribution of data */
    328:   DMSetUp(networkdm);
    329:   DMNetworkDistribute(&networkdm,0);
    330:   DMNetworkAssembleGraphStructures(networkdm);

    332:   /* We don't use these data structures anymore since they have been copied to networkdm */
    333:   if (!rank) {
    334:     PetscFree(edgelist);
    335:     PetscFree2(node,branch);
    336:   }

    338:   /* Create vectors and matrix */
    339:   DMCreateGlobalVector(networkdm,&x);
    340:   VecDuplicate(x,&b);
    341:   DMCreateMatrix(networkdm,&A);

    343:   PetscLogStagePop();

    345:   PetscLogStagePush(stage[2]);
    346:   /* Assembly system of equations */
    347:   FormOperator(networkdm,A,b);

    349:   /* Solve linear system: A x = b */
    350:   KSPCreate(PETSC_COMM_WORLD, &ksp);
    351:   KSPSetOperators(ksp, A, A);
    352:   KSPSetFromOptions(ksp);
    353:   KSPSolve(ksp, b, x);

    355:   PetscLogStagePop();

    357:   /* Free work space */
    358:   VecDestroy(&x);
    359:   VecDestroy(&b);
    360:   MatDestroy(&A);
    361:   KSPDestroy(&ksp);
    362:   DMDestroy(&networkdm);
    363:   PetscFinalize();
    364:   return ierr;
    365: }

    368: /*TEST

    370:    build:
    371:       requires: !single double define(PETSC_HAVE_ATTRIBUTEALIGNED)

    373:    test:
    374:       args: -ksp_converged_reason

    376:    test:
    377:       suffix: 2
    378:       nsize: 2
    379:       args: -petscpartitioner_type simple -pc_type asm -sub_pc_type ilu -ksp_converged_reason

    381:    test:
    382:       suffix: 3
    383:       nsize: 4
    384:       args: -petscpartitioner_type simple -pc_type asm -sub_pc_type lu -sub_pc_factor_shift_type nonzero -ksp_converged_reason

    386:    test:
    387:       suffix: graphindex
    388:       args: -n 20 -vertex_global_section_view -edge_global_section_view

    390:    test:
    391:       suffix: graphindex_2
    392:       nsize: 2
    393:       args: -petscpartitioner_type simple -n 20 -vertex_global_section_view -edge_global_section_view

    395: TEST*/