Actual source code: ex1_nest.c
1: static char help[] = "This example is based on ex1 using MATNEST. \n";
3: #include <petscdmnetwork.h>
4: #include <petscksp.h>
6: /* The topology looks like:
8: (1)
9: /|\
10: / | \
11: / | \
12: R R V
13: / |b4 \
14: b1 / (4) \ b2
15: / / \ R
16: / R R \
17: / / \ \
18: / / b5 b6 \ \
19: // \\
20: (2)--------- R -------- (3)
21: b3
23: Nodes: (1), ... (4)
24: Branches: b1, ... b6
25: Resistances: R
26: Voltage source: V
28: Additionally, there is a current source from (2) to (1).
29: */
31: /*
32: Structures containing physical data of circuit.
33: Note that no topology is defined
34: */
36: typedef struct {
37: PetscInt id; /* node id */
38: PetscScalar inj; /* current injection (A) */
39: PetscBool gr; /* grounded ? */
40: } Node;
42: typedef struct {
43: PetscInt id; /* branch id */
44: PetscScalar r; /* resistance (ohms) */
45: PetscScalar bat; /* battery (V) */
46: } Branch;
48: /*
49: read_data: this routine fills data structures with problem data.
50: This can be substituted by an external parser.
51: */
53: PetscErrorCode read_data(PetscInt *pnnode,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist)
54: {
55: PetscInt nnode, nbranch, i;
56: Branch *branch;
57: Node *node;
58: PetscInt *edgelist;
61: nnode = 4;
62: nbranch = 6;
64: PetscCalloc1(nnode,&node);
65: PetscCalloc1(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 would have that:
92: edgelist[2*i] = from node
93: edgelist[2*i + 1] = to node
94: */
96: PetscCalloc1(2*nbranch, &edgelist);
98: for (i = 0; i < nbranch; i++) {
99: switch (i) {
100: case 0:
101: edgelist[2*i] = 0;
102: edgelist[2*i + 1] = 1;
103: break;
104: case 1:
105: edgelist[2*i] = 0;
106: edgelist[2*i + 1] = 2;
107: break;
108: case 2:
109: edgelist[2*i] = 1;
110: edgelist[2*i + 1] = 2;
111: break;
112: case 3:
113: edgelist[2*i] = 0;
114: edgelist[2*i + 1] = 3;
115: break;
116: case 4:
117: edgelist[2*i] = 1;
118: edgelist[2*i + 1] = 3;
119: break;
120: case 5:
121: edgelist[2*i] = 2;
122: edgelist[2*i + 1] = 3;
123: break;
124: default:
125: break;
126: }
127: }
129: /* assign pointers */
130: *pnnode = nnode;
131: *pnbranch = nbranch;
132: *pedgelist = edgelist;
133: *pbranch = branch;
134: *pnode = node;
135: return 0;
136: }
138: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
139: {
140: Vec localb;
141: Branch *branch;
142: Node *node;
143: PetscInt e;
144: PetscInt v,vStart,vEnd;
145: PetscInt eStart, eEnd;
146: PetscBool ghost;
147: const PetscInt *cone;
148: PetscScalar *barr;
149: PetscInt lofst, lofst_to, lofst_fr;
150: PetscInt key;
151: PetscInt row[2],col[6];
152: PetscScalar val[6];
153: Mat e11, c12, c21, v22;
154: Mat **mats;
156: DMGetLocalVector(networkdm,&localb);
157: VecSet(b,0.0);
158: VecSet(localb,0.0);
160: VecGetArray(localb,&barr);
162: /* Get nested submatrices */
163: MatNestGetSubMats(A,NULL,NULL,&mats);
164: e11 = mats[0][0]; /* edges */
165: c12 = mats[0][1]; /* couplings */
166: c21 = mats[1][0]; /* couplings */
167: v22 = mats[1][1]; /* vertices */
169: /* Get vertices and edge range */
170: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
171: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
173: for (e = 0; e < eEnd; e++) {
174: DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
175: DMNetworkGetEdgeOffset(networkdm,e,&lofst);
177: DMNetworkGetConnectedVertices(networkdm,e,&cone);
178: DMNetworkGetVertexOffset(networkdm,cone[0],&lofst_fr);
179: DMNetworkGetVertexOffset(networkdm,cone[1],&lofst_to);
181: /* These are edge-edge and go to e11 */
182: row[0] = lofst;
183: col[0] = lofst; val[0] = 1;
184: MatSetValuesLocal(e11,1,row,1,col,val,INSERT_VALUES);
186: /* These are edge-vertex and go to c12 */
187: col[0] = lofst_to; val[0] = 1;
188: col[1] = lofst_fr; val[1] = -1;
189: MatSetValuesLocal(c12,1,row,2,col,val,INSERT_VALUES);
191: /* These are edge-vertex and go to c12 */
192: /* from node */
193: DMNetworkGetComponent(networkdm,cone[0],0,&key,(void**)&node,NULL);
195: if (!node->gr) {
196: row[0] = lofst_fr;
197: col[0] = lofst; val[0] = 1;
198: MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
199: }
201: /* to node */
202: DMNetworkGetComponent(networkdm,cone[1],0,&key,(void**)&node,NULL);
204: if (!node->gr) {
205: row[0] = lofst_to;
206: col[0] = lofst; val[0] = -1;
207: MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
208: }
210: /* TODO: this is not a nested vector. Need to implement nested vector */
211: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);
212: barr[lofst] = branch->bat;
213: }
215: for (v = vStart; v < vEnd; v++) {
216: DMNetworkIsGhostVertex(networkdm,v,&ghost);
217: if (!ghost) {
218: DMNetworkGetComponent(networkdm,v,0,&key,(void**)&node,NULL);
219: DMNetworkGetVertexOffset(networkdm,v,&lofst);
221: if (node->gr) {
222: row[0] = lofst;
223: col[0] = lofst; val[0] = 1;
224: MatSetValuesLocal(v22,1,row,1,col,val,INSERT_VALUES);
225: } else {
226: /* TODO: this is not a nested vector. Need to implement nested vector */
227: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);
228: barr[lofst] -= node->inj;
229: }
230: }
231: }
233: VecRestoreArray(localb,&barr);
235: DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
236: DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
237: DMRestoreLocalVector(networkdm,&localb);
239: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
240: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
241: return 0;
242: }
244: int main(int argc,char ** argv)
245: {
246: PetscInt i, nnode = 0, nbranch = 0;
247: PetscInt eStart, eEnd, vStart, vEnd;
248: PetscMPIInt size, rank;
249: DM networkdm;
250: Vec x, b;
251: Mat A;
252: KSP ksp;
253: PetscInt *edgelist = NULL;
254: PetscInt componentkey[2];
255: Node *node;
256: Branch *branch;
258: PetscInitialize(&argc,&argv,(char*)0,help);
259: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
260: MPI_Comm_size(PETSC_COMM_WORLD,&size);
262: /* "read" data only for processor 0 */
263: if (rank == 0) {
264: read_data(&nnode, &nbranch, &node, &branch, &edgelist);
265: }
267: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
268: DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
269: DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);
271: /* Set number of nodes/edges, add edge connectivity */
272: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
273: DMNetworkAddSubnetwork(networkdm,"",nbranch,edgelist,NULL);
275: /* Set up the network layout */
276: DMNetworkLayoutSetUp(networkdm);
278: /* Add network components (physical parameters of nodes and branches) and num of variables */
279: if (rank == 0) {
280: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
281: for (i = eStart; i < eEnd; i++) {
282: DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
283: }
285: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
286: for (i = vStart; i < vEnd; i++) {
287: DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
288: }
289: }
291: /* Network partitioning and distribution of data */
292: DMSetUp(networkdm);
293: DMNetworkDistribute(&networkdm,0);
295: DMNetworkAssembleGraphStructures(networkdm);
297: /* We don't use these data structures anymore since they have been copied to networkdm */
298: if (rank == 0) {
299: PetscFree(edgelist);
300: PetscFree(node);
301: PetscFree(branch);
302: }
304: DMCreateGlobalVector(networkdm,&x);
305: VecDuplicate(x,&b);
307: DMSetMatType(networkdm, MATNEST);
308: DMCreateMatrix(networkdm,&A);
310: /* Assembly system of equations */
311: FormOperator(networkdm,A,b);
313: KSPCreate(PETSC_COMM_WORLD, &ksp);
314: KSPSetOperators(ksp, A, A);
315: KSPSetFromOptions(ksp);
316: KSPSolve(ksp, b, x);
317: VecView(x, 0);
319: VecDestroy(&x);
320: VecDestroy(&b);
321: MatDestroy(&A);
322: KSPDestroy(&ksp);
323: DMDestroy(&networkdm);
324: PetscFinalize();
325: return 0;
326: }
328: /*TEST
330: build:
331: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
333: test:
334: args: -ksp_converged_reason
336: test:
337: suffix: 2
338: nsize: 2
339: args: -petscpartitioner_type simple -ksp_converged_reason
341: TEST*/