Actual source code: ex2.c
petsc-3.9.4 2018-09-11
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: {
45: PetscInt i, j, nedges = 0;
46: PetscInt *edgelist;
47: PetscInt nbat, ncurr, fr, to;
48: PetscReal *x, *y, value, xmax = 10.0; /* generate points in square */
49: PetscReal maxdist = 0.0, dist, alpha, beta, prob;
50: PetscRandom rnd;
51: Branch *branch;
52: Node *node;
53: Edge *head = NULL, *nnew= NULL, *aux= NULL;
56: PetscRandomCreate(PETSC_COMM_SELF,&rnd);
57: PetscRandomSetFromOptions(rnd);
59: PetscRandomSetSeed(rnd, seed);
60: PetscRandomSeed(rnd);
62: /* These parameters might be modified for experimentation */
63: nbat = (PetscInt)(0.1*nvertex);
64: ncurr = (PetscInt)(0.1*nvertex);
65: alpha = 0.6;
66: beta = 0.2;
68: PetscMalloc2(nvertex,&x,nvertex,&y);
70: PetscRandomSetInterval(rnd,0.0,xmax);
71: for (i=0; i<nvertex; i++) {
72: PetscRandomGetValueReal(rnd,&x[i]);
73: PetscRandomGetValueReal(rnd,&y[i]);
74: }
76: /* find maximum distance */
77: for (i=0; i<nvertex; i++) {
78: for (j=0; j<nvertex; j++) {
79: dist = findDistance(x[i],x[j],y[i],y[j]);
80: if (dist >= maxdist) maxdist = dist;
81: }
82: }
84: PetscRandomSetInterval(rnd,0.0,1.0);
85: for (i=0; i<nvertex; i++) {
86: for (j=0; j<nvertex; j++) {
87: if (j != i) {
88: dist = findDistance(x[i],x[j],y[i],y[j]);
89: prob = beta*PetscExpScalar(-dist/(maxdist*alpha));
90: PetscRandomGetValueReal(rnd,&value);
91: if (value <= prob) {
92: PetscMalloc1(1,&nnew);
93: if (head == NULL) {
94: head = nnew;
95: head->next = NULL;
96: head->n = nedges;
97: head->i = i;
98: head->j = j;
99: } else {
100: aux = head;
101: head = nnew;
102: head->n = nedges;
103: head->next = aux;
104: head->i = i;
105: head->j = j;
106: }
107: nedges += 1;
108: }
109: }
110: }
111: }
113: PetscMalloc1(2*nedges,&edgelist);
115: for (aux = head; aux; aux = aux->next) {
116: edgelist[(aux->n)*2] = aux->i;
117: edgelist[(aux->n)*2 + 1] = aux->j;
118: }
120: aux = head;
121: while (aux != NULL) {
122: nnew = aux;
123: aux = aux->next;
124: PetscFree(nnew);
125: }
127: PetscCalloc2(nvertex,&node,nedges,&branch);
128:
129: for (i = 0; i < nvertex; i++) {
130: node[i].id = i;
131: node[i].inj = 0;
132: node[i].gr = PETSC_FALSE;
133: }
135: for (i = 0; i < nedges; i++) {
136: branch[i].id = i;
137: branch[i].r = 1.0;
138: branch[i].bat = 0;
139: }
140:
141: /* Chose random node as ground voltage */
142: PetscRandomSetInterval(rnd,0.0,nvertex);
143: PetscRandomGetValueReal(rnd,&value);
144: node[(int)value].gr = PETSC_TRUE;
145:
146: /* Create random current and battery injectionsa */
147: for (i=0; i<ncurr; i++) {
148: PetscRandomSetInterval(rnd,0.0,nvertex);
149: PetscRandomGetValueReal(rnd,&value);
150: fr = edgelist[(int)value*2];
151: to = edgelist[(int)value*2 + 1];
152: node[fr].inj += 1.0;
153: node[to].inj -= 1.0;
154: }
156: for (i=0; i<nbat; i++) {
157: PetscRandomSetInterval(rnd,0.0,nedges);
158: PetscRandomGetValueReal(rnd,&value);
159: branch[(int)value].bat += 1.0;
160: }
162: PetscFree2(x,y);
163: PetscRandomDestroy(&rnd);
165: /* assign pointers */
166: *pnbranch = nedges;
167: *pedgelist = edgelist;
168: *pbranch = branch;
169: *pnode = node;
170: PetscFunctionReturn(ierr);
171: }
173: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
174: {
175: PetscErrorCode ierr;
176: Vec localb;
177: Branch *branch;
178: Node *node;
179: PetscInt e,v,vStart,vEnd,eStart, eEnd;
180: PetscInt lofst,lofst_to,lofst_fr,row[2],col[6];
181: PetscBool ghost;
182: const PetscInt *cone;
183: PetscScalar *barr,val[6];
186: DMGetLocalVector(networkdm,&localb);
187: VecSet(b,0.0);
188: VecSet(localb,0.0);
189: MatZeroEntries(A);
191: VecGetArray(localb,&barr);
193: /*
194: We can define the current as a "edge characteristic" and the voltage
195: and the voltage as a "vertex characteristic". With that, we can iterate
196: the list of edges and vertices, query the associated voltages and currents
197: and use them to write the Kirchoff equations.
198: */
200: /* Branch equations: i/r + uj - ui = battery */
201: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
202: for (e = 0; e < eEnd; e++) {
203: DMNetworkGetComponent(networkdm,e,0,NULL,(void**)&branch);
204: DMNetworkGetVariableOffset(networkdm,e,&lofst);
206: DMNetworkGetConnectedVertices(networkdm,e,&cone);
207: DMNetworkGetVariableOffset(networkdm,cone[0],&lofst_fr);
208: DMNetworkGetVariableOffset(networkdm,cone[1],&lofst_to);
210: barr[lofst] = branch->bat;
212: row[0] = lofst;
213: col[0] = lofst; val[0] = 1;
214: col[1] = lofst_to; val[1] = 1;
215: col[2] = lofst_fr; val[2] = -1;
216: MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);
218: /* from node */
219: DMNetworkGetComponent(networkdm,cone[0],0,NULL,(void**)&node);
221: if (!node->gr) {
222: row[0] = lofst_fr;
223: col[0] = lofst; val[0] = 1;
224: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
225: }
227: /* to node */
228: DMNetworkGetComponent(networkdm,cone[1],0,NULL,(void**)&node);
230: if (!node->gr) {
231: row[0] = lofst_to;
232: col[0] = lofst; val[0] = -1;
233: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
234: }
235: }
237: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
238: for (v = vStart; v < vEnd; v++) {
239: DMNetworkIsGhostVertex(networkdm,v,&ghost);
240: if (!ghost) {
241: DMNetworkGetComponent(networkdm,v,0,NULL,(void**)&node);
242: DMNetworkGetVariableOffset(networkdm,v,&lofst);
244: if (node->gr) {
245: row[0] = lofst;
246: col[0] = lofst; val[0] = 1;
247: MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
248: } else {
249: barr[lofst] -= node->inj;
250: }
251: }
252: }
254: VecRestoreArray(localb,&barr);
256: DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
257: DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
258: DMRestoreLocalVector(networkdm,&localb);
260: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
261: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
262: return(0);
263: }
265: int main(int argc,char ** argv)
266: {
267: PetscErrorCode ierr;
268: PetscInt i, nbranch = 0, eStart, eEnd, vStart, vEnd;
269: PetscInt seed = 0, nnode = 0;
270: PetscMPIInt size, rank;
271: DM networkdm;
272: Vec x, b;
273: Mat A;
274: KSP ksp;
275: PetscInt *edgelist = NULL;
276: PetscInt componentkey[2];
277: Node *node;
278: Branch *branch;
279: #if defined(PETSC_USE_LOG)
280: PetscLogStage stage[3];
281: #endif
283: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
284: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
285: MPI_Comm_size(PETSC_COMM_WORLD,&size);
287: PetscOptionsGetInt(NULL,NULL,"-seed",&seed,NULL);
289: PetscLogStageRegister("Network Creation", &stage[0]);
290: PetscLogStageRegister("DMNetwork data structures", &stage[1]);
291: PetscLogStageRegister("KSP", &stage[2]);
293: PetscLogStagePush(stage[0]);
294: /* "read" data only for processor 0 */
295: if (!rank) {
296: nnode = 100;
297: PetscOptionsGetInt(NULL,NULL,"-n",&nnode,NULL);
298: random_network(nnode, &nbranch, &node, &branch, &edgelist, seed);
299: }
300: PetscLogStagePop();
302: PetscLogStagePush(stage[1]);
303: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
304: DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
305: DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);
307: /* Set number of nodes/edges */
308: DMNetworkSetSizes(networkdm,1,0,&nnode,&nbranch,NULL,NULL);
309: /* Add edge connectivity */
310: DMNetworkSetEdgeList(networkdm,&edgelist,NULL);
311: /* Set up the network layout */
312: DMNetworkLayoutSetUp(networkdm);
314: /* Add network components: physical parameters of nodes and branches*/
315: if (!rank) {
316: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
317: for (i = eStart; i < eEnd; i++) {
318: DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart]);
319: DMNetworkAddNumVariables(networkdm,i,1);
320: }
322: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
323: for (i = vStart; i < vEnd; i++) {
324: DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart]);
325: /* Add number of variables */
326: DMNetworkAddNumVariables(networkdm,i,1);
327: }
328: }
330: /* Network partitioning and distribution of data */
331: DMSetUp(networkdm);
332: DMNetworkDistribute(&networkdm,0);
333: DMNetworkAssembleGraphStructures(networkdm);
335: /* We don't use these data structures anymore since they have been copied to networkdm */
336: if (!rank) {
337: PetscFree(edgelist);
338: PetscFree2(node,branch);
339: }
341: /* Create vectors and matrix */
342: DMCreateGlobalVector(networkdm,&x);
343: VecDuplicate(x,&b);
344: DMCreateMatrix(networkdm,&A);
346: PetscLogStagePop();
348: PetscLogStagePush(stage[2]);
349: /* Assembly system of equations */
350: FormOperator(networkdm,A,b);
352: /* Solve linear system: A x = b */
353: KSPCreate(PETSC_COMM_WORLD, &ksp);
354: KSPSetOperators(ksp, A, A);
355: KSPSetFromOptions(ksp);
356: KSPSolve(ksp, b, x);
358: PetscLogStagePop();
359:
360: /* Free work space */
361: VecDestroy(&x);
362: VecDestroy(&b);
363: MatDestroy(&A);
364: KSPDestroy(&ksp);
365: DMDestroy(&networkdm);
366: PetscFinalize();
367: return ierr;
368: }
371: /*TEST
373: build:
374: requires: !single double define(PETSC_HAVE_ATTRIBUTEALIGNED)
376: test:
377: args: -ksp_converged_reason
379: test:
380: suffix: 2
381: nsize: 2
382: args: -petscpartitioner_type simple -pc_type asm -sub_pc_type ilu -ksp_converged_reason
384: test:
385: suffix: 3
386: nsize: 4
387: args: -petscpartitioner_type simple -pc_type asm -sub_pc_type lu -sub_pc_factor_shift_type nonzero -ksp_converged_reason
389: test:
390: suffix: graphindex
391: args: -n 20 -vertex_global_section_view -edge_global_section_view
393: test:
394: suffix: graphindex_2
395: nsize: 2
396: args: -petscpartitioner_type simple -n 20 -vertex_global_section_view -edge_global_section_view
398: TEST*/