Actual source code: ex2.c
1: static char help[] = "This example is based on ex1.c, but generates a random network of chosen sizes and parameters. \n\
2: Usage: -n determines number of nodes. The nonnegative seed can be specified with the flag -seed, otherwise the program generates a random seed.\n\n";
4: /* T
5: Concepts: DMNetwork
6: Concepts: KSP
7: */
9: #include <petscdmnetwork.h>
10: #include <petscksp.h>
11: #include <time.h>
13: typedef struct {
14: PetscInt id; /* node id */
15: PetscScalar inj; /* current injection (A) */
16: PetscBool gr; /* grounded ? */
17: } Node;
19: typedef struct {
20: PetscInt id; /* branch id */
21: PetscScalar r; /* resistance (ohms) */
22: PetscScalar bat; /* battery (V) */
23: } Branch;
25: typedef struct Edge {
26: PetscInt n;
27: PetscInt i;
28: PetscInt j;
29: struct Edge *next;
30: } Edge;
32: PetscReal findDistance(PetscReal x1, PetscReal x2, PetscReal y1, PetscReal y2)
33: {
34: return PetscSqrtReal(PetscPowReal(x2-x1,2.0) + PetscPowReal(y2-y1,2.0));
35: }
37: /*
38: The algorithm for network formation is based on the paper:
39: Routing of Multipoint Connections, Bernard M. Waxman. 1988
40: */
42: PetscErrorCode random_network(PetscInt nvertex,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist,PetscInt seed)
43: {
44: PetscInt i, j, nedges = 0;
45: PetscInt *edgelist;
46: PetscInt nbat, ncurr, fr, to;
47: PetscReal *x, *y, value, xmax = 10.0; /* generate points in square */
48: PetscReal maxdist = 0.0, dist, alpha, beta, prob;
49: PetscRandom rnd;
50: Branch *branch;
51: Node *node;
52: Edge *head = NULL, *nnew= NULL, *aux= NULL;
55: PetscRandomCreate(PETSC_COMM_SELF,&rnd);
56: PetscRandomSetFromOptions(rnd);
58: PetscRandomSetSeed(rnd, seed);
59: PetscRandomSeed(rnd);
61: /* These parameters might be modified for experimentation */
62: nbat = (PetscInt)(0.1*nvertex);
63: ncurr = (PetscInt)(0.1*nvertex);
64: alpha = 0.6;
65: beta = 0.2;
67: PetscMalloc2(nvertex,&x,nvertex,&y);
69: PetscRandomSetInterval(rnd,0.0,xmax);
70: for (i=0; i<nvertex; i++) {
71: PetscRandomGetValueReal(rnd,&x[i]);
72: PetscRandomGetValueReal(rnd,&y[i]);
73: }
75: /* find maximum distance */
76: for (i=0; i<nvertex; i++) {
77: for (j=0; j<nvertex; j++) {
78: dist = findDistance(x[i],x[j],y[i],y[j]);
79: if (dist >= maxdist) maxdist = dist;
80: }
81: }
83: PetscRandomSetInterval(rnd,0.0,1.0);
84: for (i=0; i<nvertex; i++) {
85: for (j=0; j<nvertex; j++) {
86: if (j != i) {
87: dist = findDistance(x[i],x[j],y[i],y[j]);
88: prob = beta*PetscExpScalar(-dist/(maxdist*alpha));
89: PetscRandomGetValueReal(rnd,&value);
90: if (value <= prob) {
91: PetscMalloc1(1,&nnew);
92: if (head == NULL) {
93: head = nnew;
94: head->next = NULL;
95: head->n = nedges;
96: head->i = i;
97: head->j = j;
98: } else {
99: aux = head;
100: head = nnew;
101: head->n = nedges;
102: head->next = aux;
103: head->i = i;
104: head->j = j;
105: }
106: nedges += 1;
107: }
108: }
109: }
110: }
112: PetscMalloc1(2*nedges,&edgelist);
114: for (aux = head; aux; aux = aux->next) {
115: edgelist[(aux->n)*2] = aux->i;
116: edgelist[(aux->n)*2 + 1] = aux->j;
117: }
119: aux = head;
120: while (aux != NULL) {
121: nnew = aux;
122: aux = aux->next;
123: PetscFree(nnew);
124: }
126: PetscCalloc2(nvertex,&node,nedges,&branch);
128: for (i = 0; i < nvertex; i++) {
129: node[i].id = i;
130: node[i].inj = 0;
131: node[i].gr = PETSC_FALSE;
132: }
134: for (i = 0; i < nedges; i++) {
135: branch[i].id = i;
136: branch[i].r = 1.0;
137: branch[i].bat = 0;
138: }
140: /* Chose random node as ground voltage */
141: PetscRandomSetInterval(rnd,0.0,nvertex);
142: PetscRandomGetValueReal(rnd,&value);
143: node[(int)value].gr = PETSC_TRUE;
145: /* Create random current and battery injectionsa */
146: for (i=0; i<ncurr; i++) {
147: PetscRandomSetInterval(rnd,0.0,nvertex);
148: PetscRandomGetValueReal(rnd,&value);
149: fr = edgelist[(int)value*2];
150: to = edgelist[(int)value*2 + 1];
151: node[fr].inj += 1.0;
152: node[to].inj -= 1.0;
153: }
155: for (i=0; i<nbat; i++) {
156: PetscRandomSetInterval(rnd,0.0,nedges);
157: PetscRandomGetValueReal(rnd,&value);
158: branch[(int)value].bat += 1.0;
159: }
161: PetscFree2(x,y);
162: PetscRandomDestroy(&rnd);
164: /* assign pointers */
165: *pnbranch = nedges;
166: *pedgelist = edgelist;
167: *pbranch = branch;
168: *pnode = node;
169: return 0;
170: }
172: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
173: {
174: Vec localb;
175: Branch *branch;
176: Node *node;
177: PetscInt e,v,vStart,vEnd,eStart, eEnd;
178: PetscInt lofst,lofst_to,lofst_fr,row[2],col[6];
179: PetscBool ghost;
180: const PetscInt *cone;
181: PetscScalar *barr,val[6];
183: DMGetLocalVector(networkdm,&localb);
184: VecSet(b,0.0);
185: VecSet(localb,0.0);
186: MatZeroEntries(A);
188: VecGetArray(localb,&barr);
190: /*
191: We can define the current as a "edge characteristic" and the voltage
192: and the voltage as a "vertex characteristic". With that, we can iterate
193: the list of edges and vertices, query the associated voltages and currents
194: and use them to write the Kirchoff equations.
195: */
197: /* Branch equations: i/r + uj - ui = battery */
198: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
199: for (e = 0; e < eEnd; e++) {
200: DMNetworkGetComponent(networkdm,e,0,NULL,(void**)&branch,NULL);
201: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);
203: DMNetworkGetConnectedVertices(networkdm,e,&cone);
204: DMNetworkGetLocalVecOffset(networkdm,cone[0],ALL_COMPONENTS,&lofst_fr);
205: DMNetworkGetLocalVecOffset(networkdm,cone[1],ALL_COMPONENTS,&lofst_to);
207: barr[lofst] = branch->bat;
209: row[0] = lofst;
210: col[0] = lofst; val[0] = 1;
211: col[1] = lofst_to; val[1] = 1;
212: col[2] = lofst_fr; val[2] = -1;
213: MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);
215: /* from node */
216: DMNetworkGetComponent(networkdm,cone[0],0,NULL,(void**)&node,NULL);
218: if (!node->gr) {
219: row[0] = lofst_fr;
220: col[0] = lofst; val[0] = 1;
221: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
222: }
224: /* to node */
225: DMNetworkGetComponent(networkdm,cone[1],0,NULL,(void**)&node,NULL);
227: if (!node->gr) {
228: row[0] = lofst_to;
229: col[0] = lofst; val[0] = -1;
230: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
231: }
232: }
234: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
235: for (v = vStart; v < vEnd; v++) {
236: DMNetworkIsGhostVertex(networkdm,v,&ghost);
237: if (!ghost) {
238: DMNetworkGetComponent(networkdm,v,0,NULL,(void**)&node,NULL);
239: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);
241: if (node->gr) {
242: row[0] = lofst;
243: col[0] = lofst; val[0] = 1;
244: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
245: } else {
246: barr[lofst] -= node->inj;
247: }
248: }
249: }
251: VecRestoreArray(localb,&barr);
253: DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
254: DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
255: DMRestoreLocalVector(networkdm,&localb);
257: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
258: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
259: return 0;
260: }
262: int main(int argc,char ** argv)
263: {
264: PetscInt i, nbranch = 0, eStart, eEnd, vStart, vEnd;
265: PetscInt seed = 0, nnode = 0;
266: PetscMPIInt size, rank;
267: DM networkdm;
268: Vec x, b;
269: Mat A;
270: KSP ksp;
271: PetscInt *edgelist = NULL;
272: PetscInt componentkey[2];
273: Node *node;
274: Branch *branch;
275: #if defined(PETSC_USE_LOG)
276: PetscLogStage stage[3];
277: #endif
279: PetscInitialize(&argc,&argv,(char*)0,help);
280: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
281: MPI_Comm_size(PETSC_COMM_WORLD,&size);
283: PetscOptionsGetInt(NULL,NULL,"-seed",&seed,NULL);
285: PetscLogStageRegister("Network Creation", &stage[0]);
286: PetscLogStageRegister("DMNetwork data structures", &stage[1]);
287: PetscLogStageRegister("KSP", &stage[2]);
289: PetscLogStagePush(stage[0]);
290: /* "read" data only for processor 0 */
291: if (rank == 0) {
292: nnode = 100;
293: PetscOptionsGetInt(NULL,NULL,"-n",&nnode,NULL);
294: random_network(nnode, &nbranch, &node, &branch, &edgelist, seed);
295: }
296: PetscLogStagePop();
298: PetscLogStagePush(stage[1]);
299: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
300: DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
301: DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);
303: /* Set number of nodes/edges and edge connectivity */
304: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
305: DMNetworkAddSubnetwork(networkdm,"",nbranch,edgelist,NULL);
307: /* Set up the network layout */
308: DMNetworkLayoutSetUp(networkdm);
310: /* Add network components (physical parameters of nodes and branches) and num of variables */
311: if (rank == 0) {
312: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
313: for (i = eStart; i < eEnd; i++) {
314: DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
315: }
317: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
318: for (i = vStart; i < vEnd; i++) {
319: DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
320: }
321: }
323: /* Network partitioning and distribution of data */
324: DMSetUp(networkdm);
325: DMNetworkDistribute(&networkdm,0);
326: DMNetworkAssembleGraphStructures(networkdm);
328: /* We don't use these data structures anymore since they have been copied to networkdm */
329: if (rank == 0) {
330: PetscFree(edgelist);
331: PetscFree2(node,branch);
332: }
334: /* Create vectors and matrix */
335: DMCreateGlobalVector(networkdm,&x);
336: VecDuplicate(x,&b);
337: DMCreateMatrix(networkdm,&A);
339: PetscLogStagePop();
341: PetscLogStagePush(stage[2]);
342: /* Assembly system of equations */
343: FormOperator(networkdm,A,b);
345: /* Solve linear system: A x = b */
346: KSPCreate(PETSC_COMM_WORLD, &ksp);
347: KSPSetOperators(ksp, A, A);
348: KSPSetFromOptions(ksp);
349: KSPSolve(ksp, b, x);
351: PetscLogStagePop();
353: /* Free work space */
354: VecDestroy(&x);
355: VecDestroy(&b);
356: MatDestroy(&A);
357: KSPDestroy(&ksp);
358: DMDestroy(&networkdm);
359: PetscFinalize();
360: return 0;
361: }
363: /*TEST
365: build:
366: requires: !single double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
368: test:
369: args: -ksp_converged_reason
371: test:
372: suffix: 2
373: nsize: 2
374: args: -petscpartitioner_type simple -pc_type asm -sub_pc_type ilu -ksp_converged_reason
376: test:
377: suffix: 3
378: nsize: 4
379: args: -petscpartitioner_type simple -pc_type asm -sub_pc_type lu -sub_pc_factor_shift_type nonzero -ksp_converged_reason
381: test:
382: suffix: graphindex
383: args: -n 20 -vertex_global_section_view -edge_global_section_view
385: test:
386: suffix: graphindex_2
387: nsize: 2
388: args: -petscpartitioner_type simple -n 20 -vertex_global_section_view -edge_global_section_view
390: TEST*/