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; val[0] = 1./branch->r;
175: col[1] = lofst_to; val[1] = 1;
176: col[2] = lofst_fr; val[2] = -1;
177: MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);
179: /* set Node equation */
180: DMNetworkGetComponent(dmnetwork,cone[0],0,NULL,(void**)&node,NULL);
182: /* from node */
183: if (!node->gr) { /* not a boundary node */
184: row[0] = lofst_fr;
185: col[0] = lofst; val[0] = -1;
186: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
187: }
189: /* to node */
190: DMNetworkGetComponent(dmnetwork,cone[1],0,NULL,(void**)&node,NULL);
192: if (!node->gr) { /* not a boundary node */
193: row[0] = lofst_to;
194: col[0] = lofst; val[0] = 1;
195: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
196: }
197: }
199: /* set rhs b for Node equation */
200: DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
201: for (v = vStart; v < vEnd; v++) {
202: DMNetworkIsGhostVertex(dmnetwork,v,&ghost);
203: if (!ghost) {
204: DMNetworkGetComponent(dmnetwork,v,0,NULL,(void**)&node,NULL);
205: DMNetworkGetLocalVecOffset(dmnetwork,v,ALL_COMPONENTS,&lofst);
207: if (node->gr) { /* a boundary node */
208: row[0] = lofst;
209: col[0] = lofst; val[0] = 1;
210: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
211: } else { /* not a boundary node */
212: barr[lofst] += node->inj;
213: }
214: }
215: }
217: VecRestoreArray(b,&barr);
219: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
220: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
221: return 0;
222: }
224: int main(int argc,char ** argv)
225: {
226: PetscInt i, nnode = 0, nbranch = 0, eStart, eEnd, vStart, vEnd;
227: PetscMPIInt size, rank;
228: DM dmnetwork;
229: Vec x, b;
230: Mat A;
231: KSP ksp;
232: PetscInt *edgelist = NULL;
233: PetscInt componentkey[2];
234: Node *node;
235: Branch *branch;
236: PetscInt nE[1];
238: PetscInitialize(&argc,&argv,(char*)0,help);
239: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
240: MPI_Comm_size(PETSC_COMM_WORLD,&size);
242: /* "Read" data only for processor 0 */
243: if (rank == 0) {
244: read_data(&nnode, &nbranch, &node, &branch, &edgelist);
245: }
247: DMNetworkCreate(PETSC_COMM_WORLD,&dmnetwork);
248: DMNetworkRegisterComponent(dmnetwork,"nstr",sizeof(Node),&componentkey[0]);
249: DMNetworkRegisterComponent(dmnetwork,"bsrt",sizeof(Branch),&componentkey[1]);
251: /* Set local number of nodes/edges, add edge connectivity */
252: nE[0] = nbranch;
253: DMNetworkSetNumSubNetworks(dmnetwork,PETSC_DECIDE,1);
254: DMNetworkAddSubnetwork(dmnetwork,"",nE[0],edgelist,NULL);
256: /* Set up the network layout */
257: DMNetworkLayoutSetUp(dmnetwork);
259: /* Add network components (physical parameters of nodes and branches) and num of variables */
260: if (rank == 0) {
261: DMNetworkGetEdgeRange(dmnetwork,&eStart,&eEnd);
262: for (i = eStart; i < eEnd; i++) {
263: DMNetworkAddComponent(dmnetwork,i,componentkey[1],&branch[i-eStart],1);
264: }
266: DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
267: for (i = vStart; i < vEnd; i++) {
268: DMNetworkAddComponent(dmnetwork,i,componentkey[0],&node[i-vStart],1);
269: }
270: }
272: /* Network partitioning and distribution of data */
273: DMSetUp(dmnetwork);
274: DMNetworkDistribute(&dmnetwork,0);
276: /* We do not use these data structures anymore since they have been copied to dmnetwork */
277: if (rank == 0) {
278: PetscFree(edgelist);
279: PetscFree2(node,branch);
280: }
282: /* Create vectors and matrix */
283: DMCreateGlobalVector(dmnetwork,&x);
284: VecDuplicate(x,&b);
285: DMCreateMatrix(dmnetwork,&A);
287: /* Assembly system of equations */
288: FormOperator(dmnetwork,A,b);
290: /* Solve linear system: A x = b */
291: KSPCreate(PETSC_COMM_WORLD, &ksp);
292: KSPSetOperators(ksp, A, A);
293: KSPSetFromOptions(ksp);
294: KSPSolve(ksp, b, x);
295: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
297: /* Free work space */
298: VecDestroy(&x);
299: VecDestroy(&b);
300: MatDestroy(&A);
301: KSPDestroy(&ksp);
302: DMDestroy(&dmnetwork);
303: PetscFinalize();
304: return 0;
305: }
307: /*TEST
309: build:
310: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
312: test:
313: args: -ksp_monitor_short
315: test:
316: suffix: 2
317: nsize: 2
318: args: -petscpartitioner_type simple -ksp_converged_reason
320: TEST*/