Actual source code: ex3.c

  1: static char help[] = "This example tests subnetwork coupling. \n\
  2:               \n\n";

  4: #include <petscdmnetwork.h>

  6: typedef struct {
  7:   PetscInt id;
  8: } Comp0;

 10: typedef struct {
 11:   PetscScalar val;
 12: } Comp1;

 14: int main(int argc, char **argv)
 15: {
 16:   PetscMPIInt     size, rank;
 17:   DM              dmnetwork;
 18:   PetscInt        i, j, net, Nsubnet, nsubnet, ne, nv, nvar, v, ncomp, compkey0, compkey1, compkey, goffset, row;
 19:   PetscInt        numEdges[10], *edgelist[10], asvtx, bsvtx;
 20:   const PetscInt *vtx, *edges;
 21:   PetscBool       sharedv, ghost, distribute = PETSC_TRUE, test = PETSC_FALSE;
 22:   Vec             X;
 23:   Comp0           comp0;
 24:   Comp1           comp1;
 25:   PetscScalar     val;

 27:   PetscFunctionBeginUser;
 28:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 29:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 30:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 32:   /* Create a network of subnetworks */
 33:   nsubnet = 1;
 34:   if (size == 1) nsubnet = 2;

 36:   /* Create a dmnetwork and register components */
 37:   PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &dmnetwork));
 38:   PetscCall(DMNetworkRegisterComponent(dmnetwork, "comp0", sizeof(Comp0), &compkey0));
 39:   PetscCall(DMNetworkRegisterComponent(dmnetwork, "comp1", sizeof(Comp1), &compkey1));

 41:   /* Set componnet values - intentionally take rank-dependent value for test */
 42:   comp0.id  = rank;
 43:   comp1.val = 10.0 * rank;

 45:   /* Set number of subnetworks, numbers of vertices and edges over each subnetwork */
 46:   PetscCall(DMNetworkSetNumSubNetworks(dmnetwork, nsubnet, PETSC_DECIDE));
 47:   PetscCall(DMNetworkGetNumSubNetworks(dmnetwork, NULL, &Nsubnet));

 49:   /* Input subnetworks; when size>1, process[i] creates subnetwork[i] */
 50:   for (i = 0; i < Nsubnet; i++) numEdges[i] = 0;
 51:   for (i = 0; i < Nsubnet; i++) {
 52:     if (i == 0 && (size == 1 || (rank == i && size > 1))) {
 53:       numEdges[i] = 3;
 54:       PetscCall(PetscMalloc1(2 * numEdges[i], &edgelist[i]));
 55:       edgelist[i][0] = 0;
 56:       edgelist[i][1] = 2;
 57:       edgelist[i][2] = 2;
 58:       edgelist[i][3] = 1;
 59:       edgelist[i][4] = 1;
 60:       edgelist[i][5] = 3;

 62:     } else if (i == 1 && (size == 1 || (rank == i && size > 1))) {
 63:       numEdges[i] = 3;
 64:       PetscCall(PetscMalloc1(2 * numEdges[i], &edgelist[i]));
 65:       edgelist[i][0] = 0;
 66:       edgelist[i][1] = 3;
 67:       edgelist[i][2] = 3;
 68:       edgelist[i][3] = 2;
 69:       edgelist[i][4] = 2;
 70:       edgelist[i][5] = 1;

 72:     } else if (i > 1 && (size == 1 || (rank == i && size > 1))) {
 73:       numEdges[i] = 3;
 74:       PetscCall(PetscMalloc1(2 * numEdges[i], &edgelist[i]));
 75:       for (j = 0; j < numEdges[i]; j++) {
 76:         edgelist[i][2 * j]     = j;
 77:         edgelist[i][2 * j + 1] = j + 1;
 78:       }
 79:     }
 80:   }

 82:   /* Add subnetworks */
 83:   for (i = 0; i < Nsubnet; i++) {
 84:     PetscInt netNum = -1;
 85:     PetscCall(DMNetworkAddSubnetwork(dmnetwork, NULL, numEdges[i], edgelist[i], &netNum));
 86:   }

 88:   /* Add shared vertices -- all processes hold this info at current implementation */
 89:   asvtx = bsvtx = 0;
 90:   for (j = 1; j < Nsubnet; j++) {
 91:     /* vertex subnet[0].0 shares with vertex subnet[j].0 */
 92:     PetscCall(DMNetworkAddSharedVertices(dmnetwork, 0, j, 1, &asvtx, &bsvtx));
 93:   }

 95:   /* Setup the network layout */
 96:   PetscCall(DMNetworkLayoutSetUp(dmnetwork));

 98:   /* Get Subnetwork(); Add nvar=1 to subnet[0] and nvar=2 to other subnets */
 99:   for (net = 0; net < Nsubnet; net++) {
100:     PetscCall(DMNetworkGetSubnetwork(dmnetwork, net, &nv, &ne, &vtx, &edges));
101:     for (v = 0; v < nv; v++) {
102:       PetscCall(DMNetworkIsSharedVertex(dmnetwork, vtx[v], &sharedv));
103:       if (sharedv) continue;

105:       if (!net) {
106:         PetscCall(DMNetworkAddComponent(dmnetwork, vtx[v], compkey0, &comp0, 1));
107:       } else {
108:         PetscCall(DMNetworkAddComponent(dmnetwork, vtx[v], compkey1, &comp1, 2));
109:       }
110:     }
111:   }

113:   /* Add components and nvar to shared vertex -- owning and all ghost ranks must call DMNetworkAddComponent() */
114:   PetscCall(DMNetworkGetSharedVertices(dmnetwork, &nv, &vtx));
115:   for (v = 0; v < nv; v++) {
116:     PetscCall(DMNetworkAddComponent(dmnetwork, vtx[v], compkey0, &comp0, 1));
117:     PetscCall(DMNetworkAddComponent(dmnetwork, vtx[v], compkey1, &comp1, 2));
118:   }

120:   /* Enable runtime option of graph partition type -- must be called before DMSetUp() */
121:   if (size > 1) {
122:     DM               plexdm;
123:     PetscPartitioner part;
124:     PetscCall(DMNetworkGetPlex(dmnetwork, &plexdm));
125:     PetscCall(DMPlexGetPartitioner(plexdm, &part));
126:     PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE));
127:     PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_csr_alg", "mat")); /* for parmetis */
128:   }

130:   /* Setup dmnetwork */
131:   PetscCall(DMSetUp(dmnetwork));

133:   /* Redistribute the network layout; use '-distribute false' to skip */
134:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-distribute", &distribute, NULL));
135:   if (distribute) {
136:     PetscCall(DMNetworkDistribute(&dmnetwork, 0));
137:     PetscCall(DMView(dmnetwork, PETSC_VIEWER_STDOUT_WORLD));
138:   }

140:   /* Create a global vector */
141:   PetscCall(DMCreateGlobalVector(dmnetwork, &X));
142:   PetscCall(VecSet(X, 0.0));

144:   /* Set X values at shared vertex */
145:   PetscCall(DMNetworkGetSharedVertices(dmnetwork, &nv, &vtx));
146:   for (v = 0; v < nv; v++) {
147:     PetscCall(DMNetworkIsGhostVertex(dmnetwork, vtx[v], &ghost));
148:     if (ghost) continue;

150:     /* only one process holds a non-ghost vertex */
151:     PetscCall(DMNetworkGetComponent(dmnetwork, vtx[v], ALL_COMPONENTS, NULL, NULL, &nvar));
152:     PetscCall(DMNetworkGetNumComponents(dmnetwork, vtx[v], &ncomp));
153:     /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] shared v %" PetscInt_FMT ": nvar %" PetscInt_FMT ", ncomp %" PetscInt_FMT "\n",rank,vtx[v],nvar,ncomp)); */
154:     for (j = 0; j < ncomp; j++) {
155:       PetscCall(DMNetworkGetComponent(dmnetwork, vtx[v], j, &compkey, NULL, &nvar));
156:       PetscCall(DMNetworkGetGlobalVecOffset(dmnetwork, vtx[v], j, &goffset));
157:       for (i = 0; i < nvar; i++) {
158:         row = goffset + i;
159:         val = compkey + 1.0;
160:         PetscCall(VecSetValues(X, 1, &row, &val, INSERT_VALUES));
161:       }
162:     }
163:   }
164:   PetscCall(VecAssemblyBegin(X));
165:   PetscCall(VecAssemblyEnd(X));
166:   PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));

168:   /* Test DMNetworkGetSubnetwork() */
169:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_getsubnet", &test, NULL));
170:   if (test) {
171:     net = 0;
172:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-subnet", &net, NULL));
173:     PetscCall(DMNetworkGetSubnetwork(dmnetwork, net, &nv, &ne, &vtx, &edges));
174:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] subnet %" PetscInt_FMT ": nv %" PetscInt_FMT ", ne %" PetscInt_FMT "\n", rank, net, nv, ne));
175:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
176:     PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));

178:     for (i = 0; i < nv; i++) {
179:       PetscCall(DMNetworkIsGhostVertex(dmnetwork, vtx[i], &ghost));
180:       PetscCall(DMNetworkIsSharedVertex(dmnetwork, vtx[i], &sharedv));

182:       PetscCall(DMNetworkGetNumComponents(dmnetwork, vtx[i], &ncomp));
183:       if (sharedv || ghost) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  [%d] v %" PetscInt_FMT " is shared %d, is ghost %d, ncomp %" PetscInt_FMT "\n", rank, vtx[i], sharedv, ghost, ncomp));

185:       for (j = 0; j < ncomp; j++) {
186:         void *component;
187:         PetscCall(DMNetworkGetComponent(dmnetwork, vtx[i], j, &compkey, (void **)&component, NULL));
188:         if (compkey == 0) {
189:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  [%d] v %" PetscInt_FMT " compkey %" PetscInt_FMT ", mycomp0->id %" PetscInt_FMT "\n", rank, vtx[i], compkey, ((Comp0 *)component)->id));
190:         } else if (compkey == 1) {
191:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  [%d] v %" PetscInt_FMT " compkey %" PetscInt_FMT ", mycomp1->val %g\n", rank, vtx[i], compkey, (double)PetscRealPart(((Comp1 *)component)->val)));
192:         }
193:       }
194:     }
195:   }

197:   /* Free work space */
198:   PetscCall(VecDestroy(&X));
199:   for (i = 0; i < Nsubnet; i++) {
200:     if (size == 1 || rank == i) PetscCall(PetscFree(edgelist[i]));
201:   }

203:   PetscCall(DMDestroy(&dmnetwork));
204:   PetscCall(PetscFinalize());
205:   return 0;
206: }

208: /*TEST

210:    build:
211:       requires: !single double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

213:    test:
214:       args:

216:    test:
217:       suffix: 2
218:       nsize: 2
219:       args: -options_left no

221:    test:
222:       suffix: 3
223:       nsize: 4
224:       args: -options_left no

226: TEST*/