Actual source code: ex1_nest.c
1: static char help[] = "This example is based on ex1 using MATNEST. \n";
4: /* T
5: Concepts: DMNetwork
6: Concepts: KSP
7: */
9: #include <petscdmnetwork.h>
10: #include <petscksp.h>
12: /* The topology looks like:
14: (1)
15: /|\
16: / | \
17: / | \
18: R R V
19: / |b4 \
20: b1 / (4) \ b2
21: / / \ R
22: / R R \
23: / / \ \
24: / / b5 b6 \ \
25: // \\
26: (2)--------- R -------- (3)
27: b3
29: Nodes: (1), ... (4)
30: Branches: b1, ... b6
31: Resistances: R
32: Voltage source: V
34: Additionally, there is a current source from (2) to (1).
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; /* grounded ? */
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: PetscErrorCode ierr;
62: PetscInt nnode, nbranch, i;
63: Branch *branch;
64: Node *node;
65: PetscInt *edgelist;
67: nnode = 4;
68: nbranch = 6;
70: PetscCalloc1(nnode,&node);
71: PetscCalloc1(nbranch,&branch);
73: for (i = 0; i < nnode; i++) {
74: node[i].id = i;
75: node[i].inj = 0;
76: node[i].gr = PETSC_FALSE;
77: }
79: for (i = 0; i < nbranch; i++) {
80: branch[i].id = i;
81: branch[i].r = 1.0;
82: branch[i].bat = 0;
83: }
85: /*
86: Branch 1 contains a voltage source of 12.0 V
87: From node 0 to 1 there exists a current source of 4.0 A
88: Node 3 is grounded, hence the voltage is 0.
89: */
90: branch[1].bat = 12.0;
91: node[0].inj = -4.0;
92: node[1].inj = 4.0;
93: node[3].gr = PETSC_TRUE;
95: /*
96: edgelist is an array of len = 2*nbranch. that defines the
97: topology of the network. For branch i we would have that:
98: edgelist[2*i] = from node
99: edgelist[2*i + 1] = to node
100: */
102: PetscCalloc1(2*nbranch, &edgelist);
104: for (i = 0; i < nbranch; i++) {
105: switch (i) {
106: case 0:
107: edgelist[2*i] = 0;
108: edgelist[2*i + 1] = 1;
109: break;
110: case 1:
111: edgelist[2*i] = 0;
112: edgelist[2*i + 1] = 2;
113: break;
114: case 2:
115: edgelist[2*i] = 1;
116: edgelist[2*i + 1] = 2;
117: break;
118: case 3:
119: edgelist[2*i] = 0;
120: edgelist[2*i + 1] = 3;
121: break;
122: case 4:
123: edgelist[2*i] = 1;
124: edgelist[2*i + 1] = 3;
125: break;
126: case 5:
127: edgelist[2*i] = 2;
128: edgelist[2*i + 1] = 3;
129: break;
130: default:
131: break;
132: }
133: }
135: /* assign pointers */
136: *pnnode = nnode;
137: *pnbranch = nbranch;
138: *pedgelist = edgelist;
139: *pbranch = branch;
140: *pnode = node;
141: return(0);
142: }
144: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
145: {
146: PetscErrorCode ierr;
147: Vec localb;
148: Branch *branch;
149: Node *node;
150: PetscInt e;
151: PetscInt v,vStart,vEnd;
152: PetscInt eStart, eEnd;
153: PetscBool ghost;
154: const PetscInt *cone;
155: PetscScalar *barr;
156: PetscInt lofst, lofst_to, lofst_fr;
157: PetscInt key;
158: PetscInt row[2],col[6];
159: PetscScalar val[6];
160: Mat e11, c12, c21, v22;
161: Mat **mats;
164: DMGetLocalVector(networkdm,&localb);
165: VecSet(b,0.0);
166: VecSet(localb,0.0);
168: VecGetArray(localb,&barr);
170: /* Get nested submatrices */
171: MatNestGetSubMats(A,NULL,NULL,&mats);
172: e11 = mats[0][0]; /* edges */
173: c12 = mats[0][1]; /* couplings */
174: c21 = mats[1][0]; /* couplings */
175: v22 = mats[1][1]; /* vertices */
177: /* Get vertices and edge range */
178: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
179: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
181: for (e = 0; e < eEnd; e++) {
182: DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
183: DMNetworkGetEdgeOffset(networkdm,e,&lofst);
185: DMNetworkGetConnectedVertices(networkdm,e,&cone);
186: DMNetworkGetVertexOffset(networkdm,cone[0],&lofst_fr);
187: DMNetworkGetVertexOffset(networkdm,cone[1],&lofst_to);
189: /* These are edge-edge and go to e11 */
190: row[0] = lofst;
191: col[0] = lofst; val[0] = 1;
192: MatSetValuesLocal(e11,1,row,1,col,val,INSERT_VALUES);
194: /* These are edge-vertex and go to c12 */
195: col[0] = lofst_to; val[0] = 1;
196: col[1] = lofst_fr; val[1] = -1;
197: MatSetValuesLocal(c12,1,row,2,col,val,INSERT_VALUES);
199: /* These are edge-vertex and go to c12 */
200: /* from node */
201: DMNetworkGetComponent(networkdm,cone[0],0,&key,(void**)&node,NULL);
203: if (!node->gr) {
204: row[0] = lofst_fr;
205: col[0] = lofst; val[0] = 1;
206: MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
207: }
209: /* to node */
210: DMNetworkGetComponent(networkdm,cone[1],0,&key,(void**)&node,NULL);
212: if (!node->gr) {
213: row[0] = lofst_to;
214: col[0] = lofst; val[0] = -1;
215: MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
216: }
218: /* TODO: this is not a nested vector. Need to implement nested vector */
219: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);
220: barr[lofst] = branch->bat;
221: }
223: for (v = vStart; v < vEnd; v++) {
224: DMNetworkIsGhostVertex(networkdm,v,&ghost);
225: if (!ghost) {
226: DMNetworkGetComponent(networkdm,v,0,&key,(void**)&node,NULL);
227: DMNetworkGetVertexOffset(networkdm,v,&lofst);
229: if (node->gr) {
230: row[0] = lofst;
231: col[0] = lofst; val[0] = 1;
232: MatSetValuesLocal(v22,1,row,1,col,val,INSERT_VALUES);
233: } else {
234: /* TODO: this is not a nested vector. Need to implement nested vector */
235: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);
236: barr[lofst] -= node->inj;
237: }
238: }
239: }
241: VecRestoreArray(localb,&barr);
243: DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
244: DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
245: DMRestoreLocalVector(networkdm,&localb);
247: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
248: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
249: return(0);
250: }
252: int main(int argc,char ** argv)
253: {
254: PetscErrorCode ierr;
255: PetscInt i, nnode = 0, nbranch = 0;
256: PetscInt eStart, eEnd, vStart, vEnd;
257: PetscMPIInt size, rank;
258: DM networkdm;
259: Vec x, b;
260: Mat A;
261: KSP ksp;
262: PetscInt *edgelist = NULL;
263: PetscInt componentkey[2];
264: Node *node;
265: Branch *branch;
267: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
268: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
269: MPI_Comm_size(PETSC_COMM_WORLD,&size);
271: /* "read" data only for processor 0 */
272: if (!rank) {
273: read_data(&nnode, &nbranch, &node, &branch, &edgelist);
274: }
276: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
277: DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
278: DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);
280: /* Set number of nodes/edges, add edge connectivity */
281: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
282: DMNetworkAddSubnetwork(networkdm,"",nnode,nbranch,edgelist,NULL);
284: /* Set up the network layout */
285: DMNetworkLayoutSetUp(networkdm);
287: /* Add network components (physical parameters of nodes and branches) and num of variables */
288: if (!rank) {
289: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
290: for (i = eStart; i < eEnd; i++) {
291: DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
292: }
294: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
295: for (i = vStart; i < vEnd; i++) {
296: DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
297: }
298: }
300: /* Network partitioning and distribution of data */
301: DMSetUp(networkdm);
302: DMNetworkDistribute(&networkdm,0);
304: DMNetworkAssembleGraphStructures(networkdm);
306: /* We don't use these data structures anymore since they have been copied to networkdm */
307: if (!rank) {
308: PetscFree(edgelist);
309: PetscFree(node);
310: PetscFree(branch);
311: }
313: DMCreateGlobalVector(networkdm,&x);
314: VecDuplicate(x,&b);
316: DMSetMatType(networkdm, MATNEST);
317: DMCreateMatrix(networkdm,&A);
319: /* Assembly system of equations */
320: FormOperator(networkdm,A,b);
322: KSPCreate(PETSC_COMM_WORLD, &ksp);
323: KSPSetOperators(ksp, A, A);
324: KSPSetFromOptions(ksp);
325: KSPSolve(ksp, b, x);
326: VecView(x, 0);
328: VecDestroy(&x);
329: VecDestroy(&b);
330: MatDestroy(&A);
331: KSPDestroy(&ksp);
332: DMDestroy(&networkdm);
333: PetscFinalize();
334: return ierr;
335: }
338: /*TEST
340: build:
341: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
343: test:
344: args: -ksp_converged_reason
346: test:
347: suffix: 2
348: nsize: 2
349: args: -petscpartitioner_type simple -ksp_converged_reason
351: TEST*/