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: /* 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: PetscInt nnode, nbranch, i;
62: Branch *branch;
63: Node *node;
64: PetscInt *edgelist;
67: nnode = 4;
68: nbranch = 6;
70: PetscCalloc2(nnode,&node,nbranch,&branch);
72: for (i = 0; i < nnode; i++) {
73: node[i].id = i;
74: node[i].inj = 0;
75: node[i].gr = PETSC_FALSE;
76: }
78: for (i = 0; i < nbranch; i++) {
79: branch[i].id = i;
80: branch[i].r = 1.0;
81: branch[i].bat = 0;
82: }
84: /*
85: Branch 1 contains a voltage source of 12.0 V
86: From node 0 to 1 there exists a current source of 4.0 A
87: Node 3 is grounded, hence the voltage is 0.
88: */
89: branch[1].bat = 12.0;
90: node[0].inj = -4.0;
91: node[1].inj = 4.0;
92: node[3].gr = PETSC_TRUE;
94: /*
95: edgelist is an array of len = 2*nbranch that defines the
96: topology of the network. For branch i we have:
97: edgelist[2*i] = from node
98: edgelist[2*i + 1] = to node.
99: */
100: PetscCalloc1(2*nbranch, &edgelist);
102: for (i = 0; i < nbranch; i++) {
103: switch (i) {
104: case 0:
105: edgelist[2*i] = 0;
106: edgelist[2*i + 1] = 1;
107: break;
108: case 1:
109: edgelist[2*i] = 0;
110: edgelist[2*i + 1] = 2;
111: break;
112: case 2:
113: edgelist[2*i] = 1;
114: edgelist[2*i + 1] = 2;
115: break;
116: case 3:
117: edgelist[2*i] = 0;
118: edgelist[2*i + 1] = 3;
119: break;
120: case 4:
121: edgelist[2*i] = 1;
122: edgelist[2*i + 1] = 3;
123: break;
124: case 5:
125: edgelist[2*i] = 2;
126: edgelist[2*i + 1] = 3;
127: break;
128: default:
129: break;
130: }
131: }
133: /* assign pointers */
134: *pnnode = nnode;
135: *pnbranch = nbranch;
136: *pedgelist = edgelist;
137: *pbranch = branch;
138: *pnode = node;
139: return 0;
140: }
142: PetscErrorCode FormOperator(DM dmnetwork,Mat A,Vec b)
143: {
144: Branch *branch;
145: Node *node;
146: PetscInt e,v,vStart,vEnd,eStart, eEnd;
147: PetscInt lofst,lofst_to,lofst_fr,row[2],col[6];
148: PetscBool ghost;
149: const PetscInt *cone;
150: PetscScalar *barr,val[6];
152: MatZeroEntries(A);
154: VecSet(b,0.0);
155: VecGetArray(b,&barr);
157: /*
158: We define the current i as an "edge characteristic" and the voltage v as a "vertex characteristic".
159: We iterate the list of edges and vertices, query the associated voltages and currents
160: and use them to write the Kirchoff equations:
162: Branch equations: i/r + v_to - v_from = v_source (battery)
163: Node equations: sum(i_to) - sum(i_from) = i_source
164: */
165: DMNetworkGetEdgeRange(dmnetwork,&eStart,&eEnd);
166: for (e = 0; e < eEnd; e++) {
167: DMNetworkGetComponent(dmnetwork,e,0,NULL,(void**)&branch,NULL);
168: DMNetworkGetLocalVecOffset(dmnetwork,e,ALL_COMPONENTS,&lofst);
170: DMNetworkGetConnectedVertices(dmnetwork,e,&cone);
171: DMNetworkGetLocalVecOffset(dmnetwork,cone[0],ALL_COMPONENTS,&lofst_fr);
172: DMNetworkGetLocalVecOffset(dmnetwork,cone[1],ALL_COMPONENTS,&lofst_to);
174: /* set rhs b for Branch equation */
175: barr[lofst] = branch->bat;
177: /* set Branch equation */
178: row[0] = lofst;
179: col[0] = lofst; val[0] = 1./branch->r;
180: col[1] = lofst_to; val[1] = 1;
181: col[2] = lofst_fr; val[2] = -1;
182: MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);
184: /* set Node equation */
185: DMNetworkGetComponent(dmnetwork,cone[0],0,NULL,(void**)&node,NULL);
187: /* from node */
188: if (!node->gr) { /* not a boundary node */
189: row[0] = lofst_fr;
190: col[0] = lofst; val[0] = -1;
191: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
192: }
194: /* to node */
195: DMNetworkGetComponent(dmnetwork,cone[1],0,NULL,(void**)&node,NULL);
197: if (!node->gr) { /* not a boundary node */
198: row[0] = lofst_to;
199: col[0] = lofst; 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; val[0] = 1;
215: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
216: } else { /* not a boundary node */
217: barr[lofst] += node->inj;
218: }
219: }
220: }
222: VecRestoreArray(b,&barr);
224: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
225: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
226: return 0;
227: }
229: int main(int argc,char ** argv)
230: {
231: PetscInt i, nnode = 0, nbranch = 0, eStart, eEnd, vStart, vEnd;
232: PetscMPIInt size, rank;
233: DM dmnetwork;
234: Vec x, b;
235: Mat A;
236: KSP ksp;
237: PetscInt *edgelist = NULL;
238: PetscInt componentkey[2];
239: Node *node;
240: Branch *branch;
241: PetscInt nE[1];
243: PetscInitialize(&argc,&argv,(char*)0,help);
244: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
245: MPI_Comm_size(PETSC_COMM_WORLD,&size);
247: /* "Read" data only for processor 0 */
248: if (rank == 0) {
249: read_data(&nnode, &nbranch, &node, &branch, &edgelist);
250: }
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++) {
268: DMNetworkAddComponent(dmnetwork,i,componentkey[1],&branch[i-eStart],1);
269: }
271: DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
272: for (i = vStart; i < vEnd; i++) {
273: DMNetworkAddComponent(dmnetwork,i,componentkey[0],&node[i-vStart],1);
274: }
275: }
277: /* Network partitioning and distribution of data */
278: DMSetUp(dmnetwork);
279: DMNetworkDistribute(&dmnetwork,0);
281: /* We do not use these data structures anymore since they have been copied to dmnetwork */
282: if (rank == 0) {
283: PetscFree(edgelist);
284: PetscFree2(node,branch);
285: }
287: /* Create vectors and matrix */
288: DMCreateGlobalVector(dmnetwork,&x);
289: VecDuplicate(x,&b);
290: DMCreateMatrix(dmnetwork,&A);
292: /* Assembly system of equations */
293: FormOperator(dmnetwork,A,b);
295: /* Solve linear system: A x = b */
296: KSPCreate(PETSC_COMM_WORLD, &ksp);
297: KSPSetOperators(ksp, A, A);
298: KSPSetFromOptions(ksp);
299: KSPSolve(ksp, b, x);
300: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
302: /* Free work space */
303: VecDestroy(&x);
304: VecDestroy(&b);
305: MatDestroy(&A);
306: KSPDestroy(&ksp);
307: DMDestroy(&dmnetwork);
308: PetscFinalize();
309: return 0;
310: }
312: /*TEST
314: build:
315: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
317: test:
318: args: -ksp_monitor_short
320: test:
321: suffix: 2
322: nsize: 2
323: args: -petscpartitioner_type simple -ksp_converged_reason
325: TEST*/