Actual source code: ex1_nest.c
1: static char help[] = "This example is based on ex1 using MATNEST. \n";
3: /* T
4: Concepts: DMNetwork
5: Concepts: KSP
6: */
8: #include <petscdmnetwork.h>
9: #include <petscksp.h>
11: /* The topology looks like:
13: (1)
14: /|\
15: / | \
16: / | \
17: R R V
18: / |b4 \
19: b1 / (4) \ b2
20: / / \ R
21: / R R \
22: / / \ \
23: / / b5 b6 \ \
24: // \\
25: (2)--------- R -------- (3)
26: b3
28: Nodes: (1), ... (4)
29: Branches: b1, ... b6
30: Resistances: R
31: Voltage source: V
33: Additionally, there is a current source from (2) to (1).
34: */
36: /*
37: Structures containing physical data of circuit.
38: Note that no topology is defined
39: */
41: typedef struct {
42: PetscInt id; /* node id */
43: PetscScalar inj; /* current injection (A) */
44: PetscBool gr; /* grounded ? */
45: } Node;
47: typedef struct {
48: PetscInt id; /* branch id */
49: PetscScalar r; /* resistance (ohms) */
50: PetscScalar bat; /* battery (V) */
51: } Branch;
53: /*
54: read_data: this routine fills data structures with problem data.
55: This can be substituted by an external parser.
56: */
58: PetscErrorCode read_data(PetscInt *pnnode,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist)
59: {
60: PetscInt nnode, nbranch, i;
61: Branch *branch;
62: Node *node;
63: PetscInt *edgelist;
66: nnode = 4;
67: nbranch = 6;
69: PetscCalloc1(nnode,&node);
70: PetscCalloc1(nbranch,&branch);
72: for (i = 0; i < nnode; i++) {
73: node[i].id = i;
74: node[i].inj = 0;
75: node[i].gr = PETSC_FALSE;
76: }
78: for (i = 0; i < nbranch; i++) {
79: branch[i].id = i;
80: branch[i].r = 1.0;
81: branch[i].bat = 0;
82: }
84: /*
85: Branch 1 contains a voltage source of 12.0 V
86: From node 0 to 1 there exists a current source of 4.0 A
87: Node 3 is grounded, hence the voltage is 0.
88: */
89: branch[1].bat = 12.0;
90: node[0].inj = -4.0;
91: node[1].inj = 4.0;
92: node[3].gr = PETSC_TRUE;
94: /*
95: edgelist is an array of len = 2*nbranch. that defines the
96: topology of the network. For branch i we would have that:
97: edgelist[2*i] = from node
98: edgelist[2*i + 1] = to node
99: */
101: PetscCalloc1(2*nbranch, &edgelist);
103: for (i = 0; i < nbranch; i++) {
104: switch (i) {
105: case 0:
106: edgelist[2*i] = 0;
107: edgelist[2*i + 1] = 1;
108: break;
109: case 1:
110: edgelist[2*i] = 0;
111: edgelist[2*i + 1] = 2;
112: break;
113: case 2:
114: edgelist[2*i] = 1;
115: edgelist[2*i + 1] = 2;
116: break;
117: case 3:
118: edgelist[2*i] = 0;
119: edgelist[2*i + 1] = 3;
120: break;
121: case 4:
122: edgelist[2*i] = 1;
123: edgelist[2*i + 1] = 3;
124: break;
125: case 5:
126: edgelist[2*i] = 2;
127: edgelist[2*i + 1] = 3;
128: break;
129: default:
130: break;
131: }
132: }
134: /* assign pointers */
135: *pnnode = nnode;
136: *pnbranch = nbranch;
137: *pedgelist = edgelist;
138: *pbranch = branch;
139: *pnode = node;
140: return 0;
141: }
143: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
144: {
145: Vec localb;
146: Branch *branch;
147: Node *node;
148: PetscInt e;
149: PetscInt v,vStart,vEnd;
150: PetscInt eStart, eEnd;
151: PetscBool ghost;
152: const PetscInt *cone;
153: PetscScalar *barr;
154: PetscInt lofst, lofst_to, lofst_fr;
155: PetscInt key;
156: PetscInt row[2],col[6];
157: PetscScalar val[6];
158: Mat e11, c12, c21, v22;
159: Mat **mats;
161: DMGetLocalVector(networkdm,&localb);
162: VecSet(b,0.0);
163: VecSet(localb,0.0);
165: VecGetArray(localb,&barr);
167: /* Get nested submatrices */
168: MatNestGetSubMats(A,NULL,NULL,&mats);
169: e11 = mats[0][0]; /* edges */
170: c12 = mats[0][1]; /* couplings */
171: c21 = mats[1][0]; /* couplings */
172: v22 = mats[1][1]; /* vertices */
174: /* Get vertices and edge range */
175: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
176: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
178: for (e = 0; e < eEnd; e++) {
179: DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
180: DMNetworkGetEdgeOffset(networkdm,e,&lofst);
182: DMNetworkGetConnectedVertices(networkdm,e,&cone);
183: DMNetworkGetVertexOffset(networkdm,cone[0],&lofst_fr);
184: DMNetworkGetVertexOffset(networkdm,cone[1],&lofst_to);
186: /* These are edge-edge and go to e11 */
187: row[0] = lofst;
188: col[0] = lofst; val[0] = 1;
189: MatSetValuesLocal(e11,1,row,1,col,val,INSERT_VALUES);
191: /* These are edge-vertex and go to c12 */
192: col[0] = lofst_to; val[0] = 1;
193: col[1] = lofst_fr; val[1] = -1;
194: MatSetValuesLocal(c12,1,row,2,col,val,INSERT_VALUES);
196: /* These are edge-vertex and go to c12 */
197: /* from node */
198: DMNetworkGetComponent(networkdm,cone[0],0,&key,(void**)&node,NULL);
200: if (!node->gr) {
201: row[0] = lofst_fr;
202: col[0] = lofst; val[0] = 1;
203: MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
204: }
206: /* to node */
207: DMNetworkGetComponent(networkdm,cone[1],0,&key,(void**)&node,NULL);
209: if (!node->gr) {
210: row[0] = lofst_to;
211: col[0] = lofst; val[0] = -1;
212: MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
213: }
215: /* TODO: this is not a nested vector. Need to implement nested vector */
216: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);
217: barr[lofst] = branch->bat;
218: }
220: for (v = vStart; v < vEnd; v++) {
221: DMNetworkIsGhostVertex(networkdm,v,&ghost);
222: if (!ghost) {
223: DMNetworkGetComponent(networkdm,v,0,&key,(void**)&node,NULL);
224: DMNetworkGetVertexOffset(networkdm,v,&lofst);
226: if (node->gr) {
227: row[0] = lofst;
228: col[0] = lofst; val[0] = 1;
229: MatSetValuesLocal(v22,1,row,1,col,val,INSERT_VALUES);
230: } else {
231: /* TODO: this is not a nested vector. Need to implement nested vector */
232: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);
233: barr[lofst] -= node->inj;
234: }
235: }
236: }
238: VecRestoreArray(localb,&barr);
240: DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
241: DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
242: DMRestoreLocalVector(networkdm,&localb);
244: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
245: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
246: return 0;
247: }
249: int main(int argc,char ** argv)
250: {
251: PetscInt i, nnode = 0, nbranch = 0;
252: PetscInt eStart, eEnd, vStart, vEnd;
253: PetscMPIInt size, rank;
254: DM networkdm;
255: Vec x, b;
256: Mat A;
257: KSP ksp;
258: PetscInt *edgelist = NULL;
259: PetscInt componentkey[2];
260: Node *node;
261: Branch *branch;
263: PetscInitialize(&argc,&argv,(char*)0,help);
264: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
265: MPI_Comm_size(PETSC_COMM_WORLD,&size);
267: /* "read" data only for processor 0 */
268: if (rank == 0) {
269: read_data(&nnode, &nbranch, &node, &branch, &edgelist);
270: }
272: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
273: DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
274: DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);
276: /* Set number of nodes/edges, add edge connectivity */
277: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
278: DMNetworkAddSubnetwork(networkdm,"",nbranch,edgelist,NULL);
280: /* Set up the network layout */
281: DMNetworkLayoutSetUp(networkdm);
283: /* Add network components (physical parameters of nodes and branches) and num of variables */
284: if (rank == 0) {
285: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
286: for (i = eStart; i < eEnd; i++) {
287: DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
288: }
290: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
291: for (i = vStart; i < vEnd; i++) {
292: DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
293: }
294: }
296: /* Network partitioning and distribution of data */
297: DMSetUp(networkdm);
298: DMNetworkDistribute(&networkdm,0);
300: DMNetworkAssembleGraphStructures(networkdm);
302: /* We don't use these data structures anymore since they have been copied to networkdm */
303: if (rank == 0) {
304: PetscFree(edgelist);
305: PetscFree(node);
306: PetscFree(branch);
307: }
309: DMCreateGlobalVector(networkdm,&x);
310: VecDuplicate(x,&b);
312: DMSetMatType(networkdm, MATNEST);
313: DMCreateMatrix(networkdm,&A);
315: /* Assembly system of equations */
316: FormOperator(networkdm,A,b);
318: KSPCreate(PETSC_COMM_WORLD, &ksp);
319: KSPSetOperators(ksp, A, A);
320: KSPSetFromOptions(ksp);
321: KSPSolve(ksp, b, x);
322: VecView(x, 0);
324: VecDestroy(&x);
325: VecDestroy(&b);
326: MatDestroy(&A);
327: KSPDestroy(&ksp);
328: DMDestroy(&networkdm);
329: PetscFinalize();
330: return 0;
331: }
333: /*TEST
335: build:
336: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
338: test:
339: args: -ksp_converged_reason
341: test:
342: suffix: 2
343: nsize: 2
344: args: -petscpartitioner_type simple -ksp_converged_reason
346: TEST*/