Actual source code: ex1.c
1: static char help[] = "This example demonstrates the use of DMNetwork interface for solving a simple electric circuit. \n\
2: The example can be found in p.150 of 'Strang, Gilbert. Computational Science and Engineering. Wellesley, MA'.\n\n";
4: #include <petscdmnetwork.h>
5: #include <petscksp.h>
7: /* The topology looks like:
9: (0)
10: /|\
11: / | \
12: / | \
13: R R V
14: / |b3 \
15: b0 / (3) \ b1
16: / / \ R
17: / R R \
18: / / \ \
19: / / b4 b5 \ \
20: // \\
21: (1)--------- R -------- (2)
22: b2
24: Nodes: (0), ... (3)
25: Branches: b0, ... b5
26: Resistances: R
27: Voltage source: V
29: Additionally, there is a current source from (1) to (0).
30: */
32: /*
33: Structures containing physical data of circuit.
34: Note that no topology is defined
35: */
37: typedef struct {
38: PetscInt id; /* node id */
39: PetscScalar inj; /* current injection (A) */
40: PetscBool gr; /* boundary node */
41: } Node;
43: typedef struct {
44: PetscInt id; /* branch id */
45: PetscScalar r; /* resistance (ohms) */
46: PetscScalar bat; /* battery (V) */
47: } Branch;
49: /*
50: read_data: this routine fills data structures with problem data.
51: This can be substituted by an external parser.
52: */
54: PetscErrorCode read_data(PetscInt *pnnode, PetscInt *pnbranch, Node **pnode, Branch **pbranch, PetscInt **pedgelist)
55: {
56: PetscInt nnode, nbranch, i;
57: Branch *branch;
58: Node *node;
59: PetscInt *edgelist;
62: nnode = 4;
63: nbranch = 6;
65: PetscCalloc2(nnode, &node, 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 have:
92: edgelist[2*i] = from node
93: edgelist[2*i + 1] = to node.
94: */
95: PetscCalloc1(2 * nbranch, &edgelist);
97: for (i = 0; i < nbranch; i++) {
98: switch (i) {
99: case 0:
100: edgelist[2 * i] = 0;
101: edgelist[2 * i + 1] = 1;
102: break;
103: case 1:
104: edgelist[2 * i] = 0;
105: edgelist[2 * i + 1] = 2;
106: break;
107: case 2:
108: edgelist[2 * i] = 1;
109: edgelist[2 * i + 1] = 2;
110: break;
111: case 3:
112: edgelist[2 * i] = 0;
113: edgelist[2 * i + 1] = 3;
114: break;
115: case 4:
116: edgelist[2 * i] = 1;
117: edgelist[2 * i + 1] = 3;
118: break;
119: case 5:
120: edgelist[2 * i] = 2;
121: edgelist[2 * i + 1] = 3;
122: break;
123: default:
124: break;
125: }
126: }
128: /* assign pointers */
129: *pnnode = nnode;
130: *pnbranch = nbranch;
131: *pedgelist = edgelist;
132: *pbranch = branch;
133: *pnode = node;
134: return 0;
135: }
137: PetscErrorCode FormOperator(DM dmnetwork, Mat A, Vec b)
138: {
139: Branch *branch;
140: Node *node;
141: PetscInt e, v, vStart, vEnd, eStart, eEnd;
142: PetscInt lofst, lofst_to, lofst_fr, row[2], col[6];
143: PetscBool ghost;
144: const PetscInt *cone;
145: PetscScalar *barr, val[6];
147: MatZeroEntries(A);
149: VecSet(b, 0.0);
150: VecGetArray(b, &barr);
152: /*
153: We define the current i as an "edge characteristic" and the voltage v as a "vertex characteristic".
154: We iterate the list of edges and vertices, query the associated voltages and currents
155: and use them to write the Kirchoff equations:
157: Branch equations: i/r + v_to - v_from = v_source (battery)
158: Node equations: sum(i_to) - sum(i_from) = i_source
159: */
160: DMNetworkGetEdgeRange(dmnetwork, &eStart, &eEnd);
161: for (e = 0; e < eEnd; e++) {
162: DMNetworkGetComponent(dmnetwork, e, 0, NULL, (void **)&branch, NULL);
163: DMNetworkGetLocalVecOffset(dmnetwork, e, ALL_COMPONENTS, &lofst);
165: DMNetworkGetConnectedVertices(dmnetwork, e, &cone);
166: DMNetworkGetLocalVecOffset(dmnetwork, cone[0], ALL_COMPONENTS, &lofst_fr);
167: DMNetworkGetLocalVecOffset(dmnetwork, cone[1], ALL_COMPONENTS, &lofst_to);
169: /* set rhs b for Branch equation */
170: barr[lofst] = branch->bat;
172: /* set Branch equation */
173: row[0] = lofst;
174: col[0] = lofst;
175: val[0] = 1. / branch->r;
176: col[1] = lofst_to;
177: val[1] = 1;
178: col[2] = lofst_fr;
179: val[2] = -1;
180: MatSetValuesLocal(A, 1, row, 3, col, val, ADD_VALUES);
182: /* set Node equation */
183: DMNetworkGetComponent(dmnetwork, cone[0], 0, NULL, (void **)&node, NULL);
185: /* from node */
186: if (!node->gr) { /* not a boundary node */
187: row[0] = lofst_fr;
188: col[0] = lofst;
189: val[0] = -1;
190: MatSetValuesLocal(A, 1, row, 1, col, val, ADD_VALUES);
191: }
193: /* to node */
194: DMNetworkGetComponent(dmnetwork, cone[1], 0, NULL, (void **)&node, NULL);
196: if (!node->gr) { /* not a boundary node */
197: row[0] = lofst_to;
198: col[0] = lofst;
199: val[0] = 1;
200: MatSetValuesLocal(A, 1, row, 1, col, val, ADD_VALUES);
201: }
202: }
204: /* set rhs b for Node equation */
205: DMNetworkGetVertexRange(dmnetwork, &vStart, &vEnd);
206: for (v = vStart; v < vEnd; v++) {
207: DMNetworkIsGhostVertex(dmnetwork, v, &ghost);
208: if (!ghost) {
209: DMNetworkGetComponent(dmnetwork, v, 0, NULL, (void **)&node, NULL);
210: DMNetworkGetLocalVecOffset(dmnetwork, v, ALL_COMPONENTS, &lofst);
212: if (node->gr) { /* a boundary node */
213: row[0] = lofst;
214: col[0] = lofst;
215: val[0] = 1;
216: MatSetValuesLocal(A, 1, row, 1, col, val, ADD_VALUES);
217: } else { /* not a boundary node */
218: barr[lofst] += node->inj;
219: }
220: }
221: }
223: VecRestoreArray(b, &barr);
225: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
226: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
227: return 0;
228: }
230: int main(int argc, char **argv)
231: {
232: PetscInt i, nnode = 0, nbranch = 0, eStart, eEnd, vStart, vEnd;
233: PetscMPIInt size, rank;
234: DM dmnetwork;
235: Vec x, b;
236: Mat A;
237: KSP ksp;
238: PetscInt *edgelist = NULL;
239: PetscInt componentkey[2];
240: Node *node;
241: Branch *branch;
242: PetscInt nE[1];
245: PetscInitialize(&argc, &argv, (char *)0, help);
246: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
247: MPI_Comm_size(PETSC_COMM_WORLD, &size);
249: /* "Read" data only for processor 0 */
250: if (rank == 0) read_data(&nnode, &nbranch, &node, &branch, &edgelist);
252: DMNetworkCreate(PETSC_COMM_WORLD, &dmnetwork);
253: DMNetworkRegisterComponent(dmnetwork, "nstr", sizeof(Node), &componentkey[0]);
254: DMNetworkRegisterComponent(dmnetwork, "bsrt", sizeof(Branch), &componentkey[1]);
256: /* Set local number of nodes/edges, add edge connectivity */
257: nE[0] = nbranch;
258: DMNetworkSetNumSubNetworks(dmnetwork, PETSC_DECIDE, 1);
259: DMNetworkAddSubnetwork(dmnetwork, "", nE[0], edgelist, NULL);
261: /* Set up the network layout */
262: DMNetworkLayoutSetUp(dmnetwork);
264: /* Add network components (physical parameters of nodes and branches) and num of variables */
265: if (rank == 0) {
266: DMNetworkGetEdgeRange(dmnetwork, &eStart, &eEnd);
267: for (i = eStart; i < eEnd; i++) DMNetworkAddComponent(dmnetwork, i, componentkey[1], &branch[i - eStart], 1);
269: DMNetworkGetVertexRange(dmnetwork, &vStart, &vEnd);
270: for (i = vStart; i < vEnd; i++) DMNetworkAddComponent(dmnetwork, i, componentkey[0], &node[i - vStart], 1);
271: }
273: /* Network partitioning and distribution of data */
274: DMSetUp(dmnetwork);
275: DMNetworkDistribute(&dmnetwork, 0);
277: /* We do not use these data structures anymore since they have been copied to dmnetwork */
278: if (rank == 0) {
279: PetscFree(edgelist);
280: PetscFree2(node, branch);
281: }
283: /* Create vectors and matrix */
284: DMCreateGlobalVector(dmnetwork, &x);
285: VecDuplicate(x, &b);
286: DMCreateMatrix(dmnetwork, &A);
288: /* Assembly system of equations */
289: FormOperator(dmnetwork, A, b);
291: /* Solve linear system: A x = b */
292: KSPCreate(PETSC_COMM_WORLD, &ksp);
293: KSPSetOperators(ksp, A, A);
294: KSPSetFromOptions(ksp);
295: KSPSolve(ksp, b, x);
296: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
298: /* Free work space */
299: VecDestroy(&x);
300: VecDestroy(&b);
301: MatDestroy(&A);
302: KSPDestroy(&ksp);
303: DMDestroy(&dmnetwork);
304: PetscFinalize();
305: return 0;
306: }
308: /*TEST
310: build:
311: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
313: test:
314: args: -ksp_monitor_short
316: test:
317: suffix: 2
318: nsize: 2
319: args: -petscpartitioner_type simple -ksp_converged_reason
321: TEST*/