Actual source code: ex1.c
petsc-3.12.5 2020-03-29
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: /* T
5: Concepts: DMNetwork
6: Concepts: KSP
7: */
9: #include <petscdmnetwork.h>
10: #include <petscksp.h>
12: /* The topology looks like:
14: (0)
15: /|\
16: / | \
17: / | \
18: R R V
19: / |b3 \
20: b0 / (3) \ b1
21: / / \ R
22: / R R \
23: / / \ \
24: / / b4 b5 \ \
25: // \\
26: (1)--------- R -------- (2)
27: b2
29: Nodes: (0), ... (3)
30: Branches: b0, ... b5
31: Resistances: R
32: Voltage source: V
34: Additionally, there is a current source from (1) to (0).
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; /* boundary node */
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;
68: nnode = 4;
69: nbranch = 6;
71: PetscCalloc2(nnode,&node,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 have:
98: edgelist[2*i] = from node
99: edgelist[2*i + 1] = to node.
100: */
101: PetscCalloc1(2*nbranch, &edgelist);
103: for (i = 0; i < nbranch; i++) {
104: switch (i) {
105: case 0:
106: edgelist[2*i] = 0;
107: edgelist[2*i + 1] = 1;
108: break;
109: case 1:
110: edgelist[2*i] = 0;
111: edgelist[2*i + 1] = 2;
112: break;
113: case 2:
114: edgelist[2*i] = 1;
115: edgelist[2*i + 1] = 2;
116: break;
117: case 3:
118: edgelist[2*i] = 0;
119: edgelist[2*i + 1] = 3;
120: break;
121: case 4:
122: edgelist[2*i] = 1;
123: edgelist[2*i + 1] = 3;
124: break;
125: case 5:
126: edgelist[2*i] = 2;
127: edgelist[2*i + 1] = 3;
128: break;
129: default:
130: break;
131: }
132: }
134: /* assign pointers */
135: *pnnode = nnode;
136: *pnbranch = nbranch;
137: *pedgelist = edgelist;
138: *pbranch = branch;
139: *pnode = node;
140: return(0);
141: }
143: PetscErrorCode FormOperator(DM dmnetwork,Mat A,Vec b)
144: {
145: PetscErrorCode ierr;
146: Branch *branch;
147: Node *node;
148: PetscInt e,v,vStart,vEnd,eStart, eEnd;
149: PetscInt lofst,lofst_to,lofst_fr,row[2],col[6];
150: PetscBool ghost;
151: const PetscInt *cone;
152: PetscScalar *barr,val[6];
155: MatZeroEntries(A);
157: VecSet(b,0.0);
158: VecGetArray(b,&barr);
160: /*
161: We define the current i as an "edge characteristic" and the voltage v as a "vertex characteristic".
162: We iterate the list of edges and vertices, query the associated voltages and currents
163: and use them to write the Kirchoff equations:
165: Branch equations: i/r + v_to - v_from = v_source (battery)
166: Node equations: sum(i_to) - sum(i_from) = i_source
167: */
168: DMNetworkGetEdgeRange(dmnetwork,&eStart,&eEnd);
169: for (e = 0; e < eEnd; e++) {
170: DMNetworkGetComponent(dmnetwork,e,0,NULL,(void**)&branch);
171: DMNetworkGetVariableOffset(dmnetwork,e,&lofst);
173: DMNetworkGetConnectedVertices(dmnetwork,e,&cone);
174: DMNetworkGetVariableOffset(dmnetwork,cone[0],&lofst_fr);
175: DMNetworkGetVariableOffset(dmnetwork,cone[1],&lofst_to);
177: /* set rhs b for Branch equation */
178: barr[lofst] = branch->bat;
180: /* set Branch equation */
181: row[0] = lofst;
182: col[0] = lofst; val[0] = 1./branch->r;
183: col[1] = lofst_to; val[1] = 1;
184: col[2] = lofst_fr; val[2] = -1;
185: MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);
187: /* set Node equation */
188: DMNetworkGetComponent(dmnetwork,cone[0],0,NULL,(void**)&node);
190: /* from node */
191: if (!node->gr) { /* not a boundary node */
192: row[0] = lofst_fr;
193: col[0] = lofst; val[0] = -1;
194: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
195: }
197: /* to node */
198: DMNetworkGetComponent(dmnetwork,cone[1],0,NULL,(void**)&node);
200: if (!node->gr) { /* not a boundary node */
201: row[0] = lofst_to;
202: col[0] = lofst; val[0] = 1;
203: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
204: }
205: }
207: /* set rhs b for Node equation */
208: DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
209: for (v = vStart; v < vEnd; v++) {
210: DMNetworkIsGhostVertex(dmnetwork,v,&ghost);
211: if (!ghost) {
212: DMNetworkGetComponent(dmnetwork,v,0,NULL,(void**)&node);
213: DMNetworkGetVariableOffset(dmnetwork,v,&lofst);
215: if (node->gr) { /* a boundary node */
216: row[0] = lofst;
217: col[0] = lofst; val[0] = 1;
218: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
219: } else { /* not a boundary node */
220: barr[lofst] += node->inj;
221: }
222: }
223: }
225: VecRestoreArray(b,&barr);
227: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
228: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
229: return(0);
230: }
232: int main(int argc,char ** argv)
233: {
234: PetscErrorCode ierr;
235: PetscInt i, nnode = 0, nbranch = 0, eStart, eEnd, vStart, vEnd;
236: PetscMPIInt size, rank;
237: DM dmnetwork;
238: Vec x, b;
239: Mat A;
240: KSP ksp;
241: PetscInt *edgelist = NULL;
242: PetscInt componentkey[2];
243: Node *node;
244: Branch *branch;
245: PetscInt nV[1],nE[1],*edgelists[1];
247: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
248: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
249: MPI_Comm_size(PETSC_COMM_WORLD,&size);
251: /* "Read" data only for processor 0 */
252: if (!rank) {
253: read_data(&nnode, &nbranch, &node, &branch, &edgelist);
254: }
256: DMNetworkCreate(PETSC_COMM_WORLD,&dmnetwork);
257: DMNetworkRegisterComponent(dmnetwork,"nstr",sizeof(Node),&componentkey[0]);
258: DMNetworkRegisterComponent(dmnetwork,"bsrt",sizeof(Branch),&componentkey[1]);
260: /* Set local number of nodes/edges */
261: nV[0] = nnode; nE[0] = nbranch;
262: DMNetworkSetSizes(dmnetwork,1,nV,nE,0,NULL);
263: /* Add edge connectivity */
264: edgelists[0] = edgelist;
265: DMNetworkSetEdgeList(dmnetwork,edgelists,NULL);
266: /* Set up the network layout */
267: DMNetworkLayoutSetUp(dmnetwork);
269: /* Add network components: physical parameters of nodes and branches*/
270: if (!rank) {
271: DMNetworkGetEdgeRange(dmnetwork,&eStart,&eEnd);
272: for (i = eStart; i < eEnd; i++) {
273: DMNetworkAddComponent(dmnetwork,i,componentkey[1],&branch[i-eStart]);
274: /* Add number of variables */
275: DMNetworkAddNumVariables(dmnetwork,i,1);
276: }
278: DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
279: for (i = vStart; i < vEnd; i++) {
280: DMNetworkAddComponent(dmnetwork,i,componentkey[0],&node[i-vStart]);
281: /* Add number of variables */
282: DMNetworkAddNumVariables(dmnetwork,i,1);
283: }
284: }
286: /* Network partitioning and distribution of data */
287: DMSetUp(dmnetwork);
288: DMNetworkDistribute(&dmnetwork,0);
290: /* We do not use these data structures anymore since they have been copied to dmnetwork */
291: if (!rank) {
292: PetscFree(edgelist);
293: PetscFree2(node,branch);
294: }
296: /* Create vectors and matrix */
297: DMCreateGlobalVector(dmnetwork,&x);
298: VecDuplicate(x,&b);
299: DMCreateMatrix(dmnetwork,&A);
301: /* Assembly system of equations */
302: FormOperator(dmnetwork,A,b);
304: /* Solve linear system: A x = b */
305: KSPCreate(PETSC_COMM_WORLD, &ksp);
306: KSPSetOperators(ksp, A, A);
307: KSPSetFromOptions(ksp);
308: KSPSolve(ksp, b, x);
309: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
311: /* Free work space */
312: VecDestroy(&x);
313: VecDestroy(&b);
314: MatDestroy(&A);
315: KSPDestroy(&ksp);
316: DMDestroy(&dmnetwork);
317: PetscFinalize();
318: return ierr;
319: }
322: /*TEST
324: build:
325: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
327: test:
328: args: -ksp_monitor_short
330: test:
331: suffix: 2
332: nsize: 2
333: args: -petscpartitioner_type simple -ksp_converged_reason
335: TEST*/