Actual source code: ex1_nest.c
petsc-3.11.4 2019-09-28
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);
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);
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);
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: DMNetworkGetVariableOffset(networkdm,e,&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);
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: DMNetworkGetVariableOffset(networkdm,v,&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;
266: PetscInt nV[1],nE[1],*edgelists[1];
268: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
269: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
270: MPI_Comm_size(PETSC_COMM_WORLD,&size);
272: /* "read" data only for processor 0 */
273: if (!rank) {
274: read_data(&nnode, &nbranch, &node, &branch, &edgelist);
275: }
277: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
279: DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
280: DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);
283: /* Set number of nodes/edges */
284: nV[0] = nnode; nE[0] = nbranch;
285: DMNetworkSetSizes(networkdm,1,0,nV,nE,NULL,NULL);
286: /* Add edge connectivity */
287: edgelists[0] = edgelist;
288: DMNetworkSetEdgeList(networkdm,edgelists,NULL);
289: /* Set up the network layout */
290: DMNetworkLayoutSetUp(networkdm);
292: /* Add network components: physical parameters of nodes and branches*/
293: if (!rank) {
294: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
295: for (i = eStart; i < eEnd; i++) {
296: DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart]);
297: DMNetworkAddNumVariables(networkdm,i,1);
298: }
300: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
301: for (i = vStart; i < vEnd; i++) {
302: DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart]);
303: /* Add number of variables */
304: DMNetworkAddNumVariables(networkdm,i,1);
305: }
306: }
308: /* Network partitioning and distribution of data */
309: DMSetUp(networkdm);
310: DMNetworkDistribute(&networkdm,0);
312: DMNetworkAssembleGraphStructures(networkdm);
314: /* Print some info */
315: #if 0
316: PetscInt offset, goffset;
317: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
318: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
320: for (i = eStart; i < eEnd; i++) {
321: DMNetworkGetVariableOffset(networkdm,i,&offset);
322: DMNetworkGetVariableGlobalOffset(networkdm,i,&goffset);
323: PetscPrintf(PETSC_COMM_SELF,"rank[%d] edge %d - loff: %d, goff: %d .\n",rank,i,offset,goffset);
324: }
325: for (i = vStart; i < vEnd; i++) {
326: DMNetworkGetVariableOffset(networkdm,i,&offset);
327: DMNetworkGetVariableGlobalOffset(networkdm,i,&goffset);
328: PetscPrintf("rank[%d] vertex %d - loff: %d, goff: %d .\n",rank,i,offset,goffset);
329: }
330: #endif
332: /* We don't use these data structures anymore since they have been copied to networkdm */
333: if (!rank) {
334: PetscFree(edgelist);
335: PetscFree(node);
336: PetscFree(branch);
337: }
339: DMCreateGlobalVector(networkdm,&x);
340: VecDuplicate(x,&b);
342: DMSetMatType(networkdm, MATNEST);
343: DMCreateMatrix(networkdm,&A);
345: /* Assembly system of equations */
346: FormOperator(networkdm,A,b);
348: KSPCreate(PETSC_COMM_WORLD, &ksp);
349: KSPSetOperators(ksp, A, A);
350: KSPSetFromOptions(ksp);
351: KSPSolve(ksp, b, x);
352: VecView(x, 0);
354: VecDestroy(&x);
355: VecDestroy(&b);
356: MatDestroy(&A);
357: KSPDestroy(&ksp);
358: DMDestroy(&networkdm);
359: PetscFinalize();
360: return ierr;
361: }
364: /*TEST
366: build:
367: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
369: test:
370: args: -ksp_converged_reason
372: test:
373: suffix: 2
374: nsize: 2
375: args: -petscpartitioner_type simple -ksp_converged_reason
377: TEST*/