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