Actual source code: ex1_nest.c

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

  3: #include <petscdmnetwork.h>
  4: #include <petscksp.h>

  6: /* The topology looks like:

  8:             (1)
  9:             /|\
 10:            / | \
 11:           /  |  \
 12:          R   R   V
 13:         /    |b4  \
 14:     b1 /    (4)    \ b2
 15:       /    /   \    R
 16:      /   R       R   \
 17:     /  /           \  \
 18:    / / b5        b6  \ \
 19:   //                   \\
 20: (2)--------- R -------- (3)
 21:              b3

 23:   Nodes:          (1), ... (4)
 24:   Branches:       b1, ... b6
 25:   Resistances:    R
 26:   Voltage source: V

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

 31: /*
 32:   Structures containing physical data of circuit.
 33:   Note that no topology is defined
 34: */

 36: typedef struct {
 37:   PetscInt    id;  /* node id */
 38:   PetscScalar inj; /* current injection (A) */
 39:   PetscBool   gr;  /* grounded ? */
 40: } Node;

 42: typedef struct {
 43:   PetscInt    id;  /* branch id */
 44:   PetscScalar r;   /* resistance (ohms) */
 45:   PetscScalar bat; /* battery (V) */
 46: } Branch;

 48: /*
 49:   read_data: this routine fills data structures with problem data.
 50:   This can be substituted by an external parser.
 51: */

 53: PetscErrorCode read_data(PetscInt *pnnode, PetscInt *pnbranch, Node **pnode, Branch **pbranch, PetscInt **pedgelist)
 54: {
 55:   PetscInt  nnode, nbranch, i;
 56:   Branch   *branch;
 57:   Node     *node;
 58:   PetscInt *edgelist;

 61:   nnode   = 4;
 62:   nbranch = 6;

 64:   PetscCalloc1(nnode, &node);
 65:   PetscCalloc1(nbranch, &branch);

 67:   for (i = 0; i < nnode; i++) {
 68:     node[i].id  = i;
 69:     node[i].inj = 0;
 70:     node[i].gr  = PETSC_FALSE;
 71:   }

 73:   for (i = 0; i < nbranch; i++) {
 74:     branch[i].id  = i;
 75:     branch[i].r   = 1.0;
 76:     branch[i].bat = 0;
 77:   }

 79:   /*
 80:     Branch 1 contains a voltage source of 12.0 V
 81:     From node 0 to 1 there exists a current source of 4.0 A
 82:     Node 3 is grounded, hence the voltage is 0.
 83:   */
 84:   branch[1].bat = 12.0;
 85:   node[0].inj   = -4.0;
 86:   node[1].inj   = 4.0;
 87:   node[3].gr    = PETSC_TRUE;

 89:   /*
 90:     edgelist is an array of len = 2*nbranch. that defines the
 91:     topology of the network. For branch i we would have that:
 92:       edgelist[2*i]     = from node
 93:       edgelist[2*i + 1] = to node
 94:   */

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

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

129:   /* assign pointers */
130:   *pnnode    = nnode;
131:   *pnbranch  = nbranch;
132:   *pedgelist = edgelist;
133:   *pbranch   = branch;
134:   *pnode     = node;
135:   return 0;
136: }

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

156:   DMGetLocalVector(networkdm, &localb);
157:   VecSet(b, 0.0);
158:   VecSet(localb, 0.0);

160:   VecGetArray(localb, &barr);

162:   /* Get nested submatrices */
163:   MatNestGetSubMats(A, NULL, NULL, &mats);
164:   e11 = mats[0][0]; /* edges */
165:   c12 = mats[0][1]; /* couplings */
166:   c21 = mats[1][0]; /* couplings */
167:   v22 = mats[1][1]; /* vertices */

169:   /* Get vertices and edge range */
170:   DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd);
171:   DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);

173:   for (e = 0; e < eEnd; e++) {
174:     DMNetworkGetComponent(networkdm, e, 0, &key, (void **)&branch, NULL);
175:     DMNetworkGetEdgeOffset(networkdm, e, &lofst);

177:     DMNetworkGetConnectedVertices(networkdm, e, &cone);
178:     DMNetworkGetVertexOffset(networkdm, cone[0], &lofst_fr);
179:     DMNetworkGetVertexOffset(networkdm, cone[1], &lofst_to);

181:     /* These are edge-edge and go to e11 */
182:     row[0] = lofst;
183:     col[0] = lofst;
184:     val[0] = 1;
185:     MatSetValuesLocal(e11, 1, row, 1, col, val, INSERT_VALUES);

187:     /* These are edge-vertex and go to c12 */
188:     col[0] = lofst_to;
189:     val[0] = 1;
190:     col[1] = lofst_fr;
191:     val[1] = -1;
192:     MatSetValuesLocal(c12, 1, row, 2, col, val, INSERT_VALUES);

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

198:     if (!node->gr) {
199:       row[0] = lofst_fr;
200:       col[0] = lofst;
201:       val[0] = 1;
202:       MatSetValuesLocal(c21, 1, row, 1, col, val, INSERT_VALUES);
203:     }

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

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

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

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

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

239:   VecRestoreArray(localb, &barr);

241:   DMLocalToGlobalBegin(networkdm, localb, ADD_VALUES, b);
242:   DMLocalToGlobalEnd(networkdm, localb, ADD_VALUES, b);
243:   DMRestoreLocalVector(networkdm, &localb);

245:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
246:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
247:   return 0;
248: }

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

265:   PetscInitialize(&argc, &argv, (char *)0, help);
266:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
267:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

269:   /* "read" data only for processor 0 */
270:   if (rank == 0) read_data(&nnode, &nbranch, &node, &branch, &edgelist);

272:   DMNetworkCreate(PETSC_COMM_WORLD, &networkdm);
273:   DMNetworkRegisterComponent(networkdm, "nstr", sizeof(Node), &componentkey[0]);
274:   DMNetworkRegisterComponent(networkdm, "bsrt", sizeof(Branch), &componentkey[1]);

276:   /* Set number of nodes/edges, add edge connectivity */
277:   DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1);
278:   DMNetworkAddSubnetwork(networkdm, "", nbranch, edgelist, NULL);

280:   /* Set up the network layout */
281:   DMNetworkLayoutSetUp(networkdm);

283:   /* Add network components (physical parameters of nodes and branches) and num of variables */
284:   if (rank == 0) {
285:     DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd);
286:     for (i = eStart; i < eEnd; i++) DMNetworkAddComponent(networkdm, i, componentkey[1], &branch[i - eStart], 1);

288:     DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);
289:     for (i = vStart; i < vEnd; i++) DMNetworkAddComponent(networkdm, i, componentkey[0], &node[i - vStart], 1);
290:   }

292:   /* Network partitioning and distribution of data */
293:   DMSetUp(networkdm);
294:   DMNetworkDistribute(&networkdm, 0);

296:   DMNetworkAssembleGraphStructures(networkdm);

298:   /* We don't use these data structures anymore since they have been copied to networkdm */
299:   if (rank == 0) {
300:     PetscFree(edgelist);
301:     PetscFree(node);
302:     PetscFree(branch);
303:   }

305:   DMCreateGlobalVector(networkdm, &x);
306:   VecDuplicate(x, &b);

308:   DMSetMatType(networkdm, MATNEST);
309:   DMCreateMatrix(networkdm, &A);

311:   /* Assembly system of equations */
312:   FormOperator(networkdm, A, b);

314:   KSPCreate(PETSC_COMM_WORLD, &ksp);
315:   KSPSetOperators(ksp, A, A);
316:   KSPSetFromOptions(ksp);
317:   KSPSolve(ksp, b, x);
318:   VecView(x, 0);

320:   VecDestroy(&x);
321:   VecDestroy(&b);
322:   MatDestroy(&A);
323:   KSPDestroy(&ksp);
324:   DMDestroy(&networkdm);
325:   PetscFinalize();
326:   return 0;
327: }

329: /*TEST

331:    build:
332:       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

334:    test:
335:       args: -ksp_converged_reason

337:    test:
338:       suffix: 2
339:       nsize: 2
340:       args: -petscpartitioner_type simple -ksp_converged_reason

342: TEST*/