Actual source code: ex1_nest.c

      1: static char help[] = "This example is based on ex1 using MATNEST. \n";

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

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

     12: /* The topology looks like:

     14:             (1)
     15:             /|\
     16:            / | \
     17:           /  |  \
     18:          R   R   V
     19:         /    |b4  \
     20:     b1 /    (4)    \ b2
     21:       /    /   \    R
     22:      /   R       R   \
     23:     /  /           \  \
     24:    / / b5        b6  \ \
     25:   //                   \\
     26: (2)--------- R -------- (3)
     27:              b3

     29:   Nodes:          (1), ... (4)
     30:   Branches:       b1, ... b6
     31:   Resistances:    R
     32:   Voltage source: V

     34:   Additionally, there is a current source from (2) to (1).
     35: */

     37: /*
     38:   Structures containing physical data of circuit.
     39:   Note that no topology is defined
     40: */

     42: typedef struct {
     43:   PetscInt    id;  /* node id */
     44:   PetscScalar inj; /* current injection (A) */
     45:   PetscBool   gr;  /* grounded ? */
     46: } Node;

     48: typedef struct {
     49:   PetscInt    id;  /* branch id */
     50:   PetscScalar r;   /* resistance (ohms) */
     51:   PetscScalar bat; /* battery (V) */
     52: } Branch;

     54: /*
     55:   read_data: this routine fills data structures with problem data.
     56:   This can be substituted by an external parser.
     57: */

     59: PetscErrorCode read_data(PetscInt *pnnode,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist)
     60: {
     61:   PetscErrorCode    ierr;
     62:   PetscInt          nnode, nbranch, i;
     63:   Branch            *branch;
     64:   Node              *node;
     65:   PetscInt          *edgelist;

     67:   nnode   = 4;
     68:   nbranch = 6;

     70:   PetscCalloc1(nnode,&node);
     71:   PetscCalloc1(nbranch,&branch);

     73:   for (i = 0; i < nnode; i++) {
     74:     node[i].id  = i;
     75:     node[i].inj = 0;
     76:     node[i].gr = PETSC_FALSE;
     77:   }

     79:   for (i = 0; i < nbranch; i++) {
     80:     branch[i].id  = i;
     81:     branch[i].r   = 1.0;
     82:     branch[i].bat = 0;
     83:   }

     85:   /*
     86:     Branch 1 contains a voltage source of 12.0 V
     87:     From node 0 to 1 there exists a current source of 4.0 A
     88:     Node 3 is grounded, hence the voltage is 0.
     89:   */
     90:   branch[1].bat = 12.0;
     91:   node[0].inj   = -4.0;
     92:   node[1].inj   =  4.0;
     93:   node[3].gr    =  PETSC_TRUE;

     95:   /*
     96:     edgelist is an array of len = 2*nbranch. that defines the
     97:     topology of the network. For branch i we would have that:
     98:       edgelist[2*i]     = from node
     99:       edgelist[2*i + 1] = to node
    100:   */

    102:   PetscCalloc1(2*nbranch, &edgelist);

    104:   for (i = 0; i < nbranch; i++) {
    105:     switch (i) {
    106:       case 0:
    107:         edgelist[2*i]     = 0;
    108:         edgelist[2*i + 1] = 1;
    109:         break;
    110:       case 1:
    111:         edgelist[2*i]     = 0;
    112:         edgelist[2*i + 1] = 2;
    113:         break;
    114:       case 2:
    115:         edgelist[2*i]     = 1;
    116:         edgelist[2*i + 1] = 2;
    117:         break;
    118:       case 3:
    119:         edgelist[2*i]     = 0;
    120:         edgelist[2*i + 1] = 3;
    121:         break;
    122:       case 4:
    123:         edgelist[2*i]     = 1;
    124:         edgelist[2*i + 1] = 3;
    125:         break;
    126:       case 5:
    127:         edgelist[2*i]     = 2;
    128:         edgelist[2*i + 1] = 3;
    129:         break;
    130:       default:
    131:         break;
    132:     }
    133:   }

    135:   /* assign pointers */
    136:   *pnnode    = nnode;
    137:   *pnbranch  = nbranch;
    138:   *pedgelist = edgelist;
    139:   *pbranch   = branch;
    140:   *pnode     = node;
    141:   return(0);
    142: }

    144: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
    145: {
    146:   PetscErrorCode    ierr;
    147:   Vec               localb;
    148:   Branch            *branch;
    149:   Node              *node;
    150:   PetscInt          e;
    151:   PetscInt          v,vStart,vEnd;
    152:   PetscInt          eStart, eEnd;
    153:   PetscBool         ghost;
    154:   const PetscInt    *cone;
    155:   PetscScalar       *barr;
    156:   PetscInt          lofst, lofst_to, lofst_fr;
    157:   PetscInt          key;
    158:   PetscInt          row[2],col[6];
    159:   PetscScalar       val[6];
    160:   Mat               e11, c12, c21, v22;
    161:   Mat               **mats;

    164:   DMGetLocalVector(networkdm,&localb);
    165:   VecSet(b,0.0);
    166:   VecSet(localb,0.0);

    168:   VecGetArray(localb,&barr);

    170:   /* Get nested submatrices */
    171:   MatNestGetSubMats(A,NULL,NULL,&mats);
    172:   e11 = mats[0][0];  /* edges */
    173:   c12 = mats[0][1];  /* couplings */
    174:   c21 = mats[1][0];  /* couplings */
    175:   v22 = mats[1][1];  /* vertices */

    177:   /* Get vertices and edge range */
    178:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    179:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

    181:   for (e = 0; e < eEnd; e++) {
    182:     DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
    183:     DMNetworkGetEdgeOffset(networkdm,e,&lofst);

    185:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
    186:     DMNetworkGetVertexOffset(networkdm,cone[0],&lofst_fr);
    187:     DMNetworkGetVertexOffset(networkdm,cone[1],&lofst_to);

    189:     /* These are edge-edge and go to e11 */
    190:     row[0] = lofst;
    191:     col[0] = lofst;     val[0] =  1;
    192:     MatSetValuesLocal(e11,1,row,1,col,val,INSERT_VALUES);

    194:     /* These are edge-vertex and go to c12 */
    195:     col[0] = lofst_to;  val[0] =  1;
    196:     col[1] = lofst_fr;  val[1] = -1;
    197:     MatSetValuesLocal(c12,1,row,2,col,val,INSERT_VALUES);

    199:     /* These are edge-vertex and go to c12 */
    200:     /* from node */
    201:     DMNetworkGetComponent(networkdm,cone[0],0,&key,(void**)&node,NULL);

    203:     if (!node->gr) {
    204:       row[0] = lofst_fr;
    205:       col[0] = lofst;   val[0] =  1;
    206:       MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
    207:     }

    209:     /* to node */
    210:     DMNetworkGetComponent(networkdm,cone[1],0,&key,(void**)&node,NULL);

    212:     if (!node->gr) {
    213:       row[0] = lofst_to;
    214:       col[0] = lofst;   val[0] =  -1;
    215:       MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
    216:     }

    218:     /* TODO: this is not a nested vector. Need to implement nested vector */
    219:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);
    220:     barr[lofst] = branch->bat;
    221:   }

    223:   for (v = vStart; v < vEnd; v++) {
    224:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
    225:     if (!ghost) {
    226:       DMNetworkGetComponent(networkdm,v,0,&key,(void**)&node,NULL);
    227:       DMNetworkGetVertexOffset(networkdm,v,&lofst);

    229:       if (node->gr) {
    230:         row[0] = lofst;
    231:         col[0] = lofst;   val[0] =  1;
    232:         MatSetValuesLocal(v22,1,row,1,col,val,INSERT_VALUES);
    233:       } else {
    234:         /* TODO: this is not a nested vector. Need to implement nested vector */
    235:         DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);
    236:         barr[lofst] -= node->inj;
    237:       }
    238:     }
    239:   }

    241:   VecRestoreArray(localb,&barr);

    243:   DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
    244:   DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
    245:   DMRestoreLocalVector(networkdm,&localb);

    247:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    248:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    249:   return(0);
    250: }

    252: int main(int argc,char ** argv)
    253: {
    254:   PetscErrorCode    ierr;
    255:   PetscInt          i, nnode = 0, nbranch = 0;
    256:   PetscInt          eStart, eEnd, vStart, vEnd;
    257:   PetscMPIInt       size, rank;
    258:   DM                networkdm;
    259:   Vec               x, b;
    260:   Mat               A;
    261:   KSP               ksp;
    262:   PetscInt          *edgelist = NULL;
    263:   PetscInt          componentkey[2];
    264:   Node              *node;
    265:   Branch            *branch;

    267:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
    268:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    269:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

    271:   /* "read" data only for processor 0 */
    272:   if (!rank) {
    273:     read_data(&nnode, &nbranch, &node, &branch, &edgelist);
    274:   }

    276:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
    277:   DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
    278:   DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);

    280:   /* Set number of nodes/edges, add edge connectivity */
    281:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
    282:   DMNetworkAddSubnetwork(networkdm,"",nnode,nbranch,edgelist,NULL);

    284:   /* Set up the network layout */
    285:   DMNetworkLayoutSetUp(networkdm);

    287:   /* Add network components (physical parameters of nodes and branches) and num of variables */
    288:   if (!rank) {
    289:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    290:     for (i = eStart; i < eEnd; i++) {
    291:       DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
    292:     }

    294:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    295:     for (i = vStart; i < vEnd; i++) {
    296:       DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
    297:     }
    298:   }

    300:   /* Network partitioning and distribution of data */
    301:   DMSetUp(networkdm);
    302:   DMNetworkDistribute(&networkdm,0);

    304:   DMNetworkAssembleGraphStructures(networkdm);

    306:   /* We don't use these data structures anymore since they have been copied to networkdm */
    307:   if (!rank) {
    308:     PetscFree(edgelist);
    309:     PetscFree(node);
    310:     PetscFree(branch);
    311:   }

    313:   DMCreateGlobalVector(networkdm,&x);
    314:   VecDuplicate(x,&b);

    316:   DMSetMatType(networkdm, MATNEST);
    317:   DMCreateMatrix(networkdm,&A);

    319:   /* Assembly system of equations */
    320:   FormOperator(networkdm,A,b);

    322:   KSPCreate(PETSC_COMM_WORLD, &ksp);
    323:   KSPSetOperators(ksp, A, A);
    324:   KSPSetFromOptions(ksp);
    325:   KSPSolve(ksp, b, x);
    326:   VecView(x, 0);

    328:   VecDestroy(&x);
    329:   VecDestroy(&b);
    330:   MatDestroy(&A);
    331:   KSPDestroy(&ksp);
    332:   DMDestroy(&networkdm);
    333:   PetscFinalize();
    334:   return ierr;
    335: }

    338: /*TEST

    340:    build:
    341:       requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)

    343:    test:
    344:       args: -ksp_converged_reason

    346:    test:
    347:       suffix: 2
    348:       nsize: 2
    349:       args: -petscpartitioner_type simple -ksp_converged_reason

    351: TEST*/