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*/