1 static char help[] = "This example demonstrates the use of DMNetwork interface with subnetworks for solving a coupled nonlinear \n\ 2 electric power grid and water pipe problem.\n\ 3 The available solver options are in the ex1options file \n\ 4 and the data files are in the datafiles of subdirectories.\n\ 5 This example shows the use of subnetwork feature in DMNetwork. \n\ 6 Run this program: mpiexec -n <n> ./ex1 \n\\n"; 7 8 /* T 9 Concepts: DMNetwork 10 Concepts: PETSc SNES solver 11 */ 12 13 #include "power/power.h" 14 #include "water/water.h" 15 16 typedef struct{ 17 UserCtx_Power appctx_power; 18 AppCtx_Water appctx_water; 19 PetscInt subsnes_id; /* snes solver id */ 20 PetscInt it; /* iteration number */ 21 Vec localXold; /* store previous solution, used by FormFunction_Dummy() */ 22 } UserCtx; 23 24 PetscErrorCode UserMonitor(SNES snes,PetscInt its,PetscReal fnorm ,void *appctx) 25 { 26 PetscErrorCode ierr; 27 UserCtx *user = (UserCtx*)appctx; 28 Vec X,localXold = user->localXold; 29 DM networkdm; 30 PetscMPIInt rank; 31 MPI_Comm comm; 32 33 PetscFunctionBegin; 34 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 35 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 36 #if 0 37 if (rank == 0) { 38 PetscInt subsnes_id = user->subsnes_id; 39 if (subsnes_id == 2) { 40 ierr = PetscPrintf(PETSC_COMM_SELF," it %D, subsnes_id %D, fnorm %g\n",user->it,user->subsnes_id,(double)fnorm);CHKERRQ(ierr); 41 } else { 42 ierr = PetscPrintf(PETSC_COMM_SELF," subsnes_id %D, fnorm %g\n",user->subsnes_id,(double)fnorm);CHKERRQ(ierr); 43 } 44 } 45 #endif 46 ierr = SNESGetSolution(snes,&X);CHKERRQ(ierr); 47 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 48 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localXold);CHKERRQ(ierr); 49 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localXold);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 PetscErrorCode FormJacobian_subPower(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx) 54 { 55 PetscErrorCode ierr; 56 DM networkdm; 57 Vec localX; 58 PetscInt nv,ne,i,j,offset,nvar,row; 59 const PetscInt *vtx,*edges; 60 PetscBool ghostvtex; 61 PetscScalar one = 1.0; 62 PetscMPIInt rank; 63 MPI_Comm comm; 64 65 PetscFunctionBegin; 66 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 67 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 68 69 ierr = PetscObjectGetComm((PetscObject)networkdm,&comm);CHKERRQ(ierr); 70 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 71 72 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 73 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 74 75 ierr = MatZeroEntries(J);CHKERRQ(ierr); 76 77 /* Power subnetwork: copied from power/FormJacobian_Power() */ 78 ierr = DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 79 ierr = FormJacobian_Power_private(networkdm,localX,J,nv,ne,vtx,edges,appctx);CHKERRQ(ierr); 80 81 /* Water subnetwork: Identity */ 82 ierr = DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 83 for (i=0; i<nv; i++) { 84 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr); 85 if (ghostvtex) continue; 86 87 ierr = DMNetworkGetGlobalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);CHKERRQ(ierr); 88 ierr = DMNetworkGetComponent(networkdm,vtx[i],ALL_COMPONENTS,NULL,NULL,&nvar);CHKERRQ(ierr); 89 for (j=0; j<nvar; j++) { 90 row = offset + j; 91 ierr = MatSetValues(J,1,&row,1,&row,&one,ADD_VALUES);CHKERRQ(ierr); 92 } 93 } 94 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 95 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 96 97 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 /* Dummy equation localF(X) = localX - localXold */ 102 PetscErrorCode FormFunction_Dummy(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx) 103 { 104 PetscErrorCode ierr; 105 const PetscScalar *xarr,*xoldarr; 106 PetscScalar *farr; 107 PetscInt i,j,offset,nvar; 108 PetscBool ghostvtex; 109 UserCtx *user = (UserCtx*)appctx; 110 Vec localXold = user->localXold; 111 112 PetscFunctionBegin; 113 ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); 114 ierr = VecGetArrayRead(localXold,&xoldarr);CHKERRQ(ierr); 115 ierr = VecGetArray(localF,&farr);CHKERRQ(ierr); 116 117 for (i=0; i<nv; i++) { 118 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr); 119 if (ghostvtex) continue; 120 121 ierr = DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);CHKERRQ(ierr); 122 ierr = DMNetworkGetComponent(networkdm,vtx[i],ALL_COMPONENTS,NULL,NULL,&nvar);CHKERRQ(ierr); 123 for (j=0; j<nvar; j++) { 124 farr[offset+j] = xarr[offset+j] - xoldarr[offset+j]; 125 } 126 } 127 128 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 129 ierr = VecRestoreArrayRead(localXold,&xoldarr);CHKERRQ(ierr); 130 ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 134 PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *appctx) 135 { 136 PetscErrorCode ierr; 137 DM networkdm; 138 Vec localX,localF; 139 PetscInt nv,ne,v; 140 const PetscInt *vtx,*edges; 141 PetscMPIInt rank; 142 MPI_Comm comm; 143 UserCtx *user = (UserCtx*)appctx; 144 UserCtx_Power appctx_power = (*user).appctx_power; 145 AppCtx_Water appctx_water = (*user).appctx_water; 146 147 PetscFunctionBegin; 148 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 149 ierr = PetscObjectGetComm((PetscObject)networkdm,&comm);CHKERRQ(ierr); 150 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 151 152 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 153 ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr); 154 ierr = VecSet(F,0.0);CHKERRQ(ierr); 155 ierr = VecSet(localF,0.0);CHKERRQ(ierr); 156 157 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 158 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 159 160 /* Form Function for power subnetwork */ 161 ierr = DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 162 if (user->subsnes_id == 1) { /* snes_water only */ 163 ierr = FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);CHKERRQ(ierr); 164 } else { 165 ierr = FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,&appctx_power);CHKERRQ(ierr); 166 } 167 168 /* Form Function for water subnetwork */ 169 ierr = DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 170 if (user->subsnes_id == 0) { /* snes_power only */ 171 ierr = FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);CHKERRQ(ierr); 172 } else { 173 ierr = FormFunction_Water(networkdm,localX,localF,nv,ne,vtx,edges,NULL);CHKERRQ(ierr); 174 } 175 176 /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */ 177 ierr = DMNetworkGetSharedVertices(networkdm,&nv,&vtx);CHKERRQ(ierr); 178 for (v=0; v<nv; v++) { 179 PetscInt key,ncomp,nvar,nconnedges,k,e,keye,goffset[3]; 180 void* component; 181 const PetscInt *connedges; 182 183 ierr = DMNetworkGetComponent(networkdm,vtx[v],ALL_COMPONENTS,NULL,NULL,&nvar);CHKERRQ(ierr); 184 ierr = DMNetworkGetNumComponents(networkdm,vtx[v],&ncomp);CHKERRQ(ierr); 185 /* printf(" [%d] coupling vertex[%D]: v %D, ncomp %D; nvar %D\n",rank,v,vtx[v], ncomp,nvar); */ 186 187 for (k=0; k<ncomp; k++) { 188 ierr = DMNetworkGetComponent(networkdm,vtx[v],k,&key,&component,&nvar);CHKERRQ(ierr); 189 ierr = DMNetworkGetGlobalVecOffset(networkdm,vtx[v],k,&goffset[k]);CHKERRQ(ierr); 190 191 /* Verify the coupling vertex is a powernet load vertex or a water vertex */ 192 switch (k) { 193 case 0: 194 if (key != appctx_power.compkey_bus || nvar != 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"key %D not a power bus vertex or nvar %D != 2",key,nvar); 195 break; 196 case 1: 197 if (key != appctx_power.compkey_load || nvar != 0 || goffset[1] != goffset[0]+2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power load vertex"); 198 break; 199 case 2: 200 if (key != appctx_water.compkey_vtx || nvar != 1 || goffset[2] != goffset[1]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex"); 201 break; 202 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %D is wrong",k); 203 } 204 /* printf(" [%d] coupling vertex[%D]: key %D; nvar %D, goffset %D\n",rank,v,key,nvar,goffset[k]); */ 205 } 206 207 /* Get its supporting edges */ 208 ierr = DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);CHKERRQ(ierr); 209 /* printf("\n[%d] coupling vertex: nconnedges %D\n",rank,nconnedges);CHKERRQ(ierr); */ 210 for (k=0; k<nconnedges; k++) { 211 e = connedges[k]; 212 ierr = DMNetworkGetNumComponents(networkdm,e,&ncomp);CHKERRQ(ierr); 213 /* printf("\n [%d] connected edge[%D]=%D has ncomp %D\n",rank,k,e,ncomp); */ 214 ierr = DMNetworkGetComponent(networkdm,e,0,&keye,&component,NULL);CHKERRQ(ierr); 215 if (keye == appctx_water.compkey_edge) { /* water_compkey_edge */ 216 EDGE_Water edge=(EDGE_Water)component; 217 if (edge->type == EDGE_TYPE_PUMP) { 218 /* printf(" connected edge[%D]=%D has keye=%D, is appctx_water.compkey_edge with EDGE_TYPE_PUMP\n",k,e,keye); */ 219 } 220 } else { /* ower->compkey_branch */ 221 if (keye != appctx_power.compkey_branch) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power branch"); 222 } 223 } 224 } 225 226 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 227 228 ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 229 ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 230 ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr); 231 #if 0 232 if (rank == 0) printf("F:\n"); 233 ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 234 #endif 235 PetscFunctionReturn(0); 236 } 237 238 PetscErrorCode SetInitialGuess(DM networkdm,Vec X,void* appctx) 239 { 240 PetscErrorCode ierr; 241 PetscInt nv,ne,i,j,ncomp,offset,key; 242 const PetscInt *vtx,*edges; 243 UserCtx *user = (UserCtx*)appctx; 244 Vec localX = user->localXold; 245 UserCtx_Power appctx_power = (*user).appctx_power; 246 AppCtx_Water appctx_water = (*user).appctx_water; 247 PetscBool ghost; 248 PetscScalar *xarr; 249 VERTEX_Power bus; 250 VERTEX_Water vertex; 251 void* component; 252 GEN gen; 253 254 PetscFunctionBegin; 255 ierr = VecSet(X,0.0);CHKERRQ(ierr); 256 ierr = VecSet(localX,0.0);CHKERRQ(ierr); 257 258 /* Set initial guess for power subnetwork */ 259 ierr = DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 260 ierr = SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,&appctx_power);CHKERRQ(ierr); 261 262 /* Set initial guess for water subnetwork */ 263 ierr = DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 264 ierr = SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);CHKERRQ(ierr); 265 266 /* Set initial guess at the coupling vertex */ 267 ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr); 268 ierr = DMNetworkGetSharedVertices(networkdm,&nv,&vtx);CHKERRQ(ierr); 269 for (i=0; i<nv; i++) { 270 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);CHKERRQ(ierr); 271 if (ghost) continue; 272 273 ierr = DMNetworkGetNumComponents(networkdm,vtx[i],&ncomp);CHKERRQ(ierr); 274 for (j=0; j<ncomp; j++) { 275 ierr = DMNetworkGetLocalVecOffset(networkdm,vtx[i],j,&offset);CHKERRQ(ierr); 276 ierr = DMNetworkGetComponent(networkdm,vtx[i],j,&key,(void**)&component,NULL);CHKERRQ(ierr); 277 if (key == appctx_power.compkey_bus) { 278 bus = (VERTEX_Power)(component); 279 xarr[offset] = bus->va*PETSC_PI/180.0; 280 xarr[offset+1] = bus->vm; 281 } else if (key == appctx_power.compkey_gen) { 282 gen = (GEN)(component); 283 if (!gen->status) continue; 284 xarr[offset+1] = gen->vs; 285 } else if (key == appctx_water.compkey_vtx) { 286 vertex = (VERTEX_Water)(component); 287 if (vertex->type == VERTEX_TYPE_JUNCTION) { 288 xarr[offset] = 100; 289 } else if (vertex->type == VERTEX_TYPE_RESERVOIR) { 290 xarr[offset] = vertex->res.head; 291 } else { 292 xarr[offset] = vertex->tank.initlvl + vertex->tank.elev; 293 } 294 } 295 } 296 } 297 ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr); 298 299 ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 300 ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 int main(int argc,char **argv) 305 { 306 PetscErrorCode ierr; 307 DM networkdm; 308 PetscLogStage stage[4]; 309 PetscMPIInt rank,size; 310 PetscInt Nsubnet=2,numVertices[2],numEdges[2],i,j,nv,ne,it_max=10; 311 const PetscInt *vtx,*edges; 312 Vec X,F; 313 SNES snes,snes_power,snes_water; 314 Mat Jac; 315 PetscBool ghost,viewJ=PETSC_FALSE,viewX=PETSC_FALSE,viewDM=PETSC_FALSE,test=PETSC_FALSE,distribute=PETSC_TRUE,flg; 316 UserCtx user; 317 SNESConvergedReason reason; 318 319 /* Power subnetwork */ 320 UserCtx_Power *appctx_power = &user.appctx_power; 321 char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m"; 322 PFDATA *pfdata = NULL; 323 PetscInt genj,loadj,*edgelist_power = NULL,power_netnum; 324 PetscScalar Sbase = 0.0; 325 326 /* Water subnetwork */ 327 AppCtx_Water *appctx_water = &user.appctx_water; 328 WATERDATA *waterdata = NULL; 329 char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp"; 330 PetscInt *edgelist_water = NULL,water_netnum; 331 332 /* Shared vertices between subnetworks */ 333 PetscInt power_svtx,water_svtx; 334 335 ierr = PetscInitialize(&argc,&argv,"ex1options",help);if (ierr) return ierr; 336 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 337 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 338 339 /* (1) Read Data - Only rank 0 reads the data */ 340 /*--------------------------------------------*/ 341 ierr = PetscLogStageRegister("Read Data",&stage[0]);CHKERRQ(ierr); 342 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 343 344 for (i=0; i<Nsubnet; i++) { 345 numVertices[i] = 0; 346 numEdges[i] = 0; 347 } 348 349 /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */ 350 /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */ 351 ierr = PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);CHKERRQ(ierr); 352 ierr = PetscNew(&pfdata);CHKERRQ(ierr); 353 ierr = PFReadMatPowerData(pfdata,pfdata_file);CHKERRQ(ierr); 354 if (rank == 0) { 355 ierr = PetscPrintf(PETSC_COMM_SELF,"Power network: nb = %D, ngen = %D, nload = %D, nbranch = %D\n",pfdata->nbus,pfdata->ngen,pfdata->nload,pfdata->nbranch);CHKERRQ(ierr); 356 } 357 Sbase = pfdata->sbase; 358 if (rank == 0) { /* proc[0] will create Electric Power Grid */ 359 numEdges[0] = pfdata->nbranch; 360 numVertices[0] = pfdata->nbus; 361 362 ierr = PetscMalloc1(2*numEdges[0],&edgelist_power);CHKERRQ(ierr); 363 ierr = GetListofEdges_Power(pfdata,edgelist_power);CHKERRQ(ierr); 364 } 365 /* Broadcast power Sbase to all processors */ 366 ierr = MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRMPI(ierr); 367 appctx_power->Sbase = Sbase; 368 appctx_power->jac_error = PETSC_FALSE; 369 /* If external option activated. Introduce error in jacobian */ 370 ierr = PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);CHKERRQ(ierr); 371 372 /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */ 373 /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */ 374 ierr = PetscNew(&waterdata);CHKERRQ(ierr); 375 ierr = PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,PETSC_MAX_PATH_LEN-1,NULL);CHKERRQ(ierr); 376 ierr = WaterReadData(waterdata,waterdata_file);CHKERRQ(ierr); 377 if (size == 1 || (size > 1 && rank == 1)) { 378 ierr = PetscCalloc1(2*waterdata->nedge,&edgelist_water);CHKERRQ(ierr); 379 ierr = GetListofEdges_Water(waterdata,edgelist_water);CHKERRQ(ierr); 380 numEdges[1] = waterdata->nedge; 381 numVertices[1] = waterdata->nvertex; 382 } 383 PetscLogStagePop(); 384 385 /* (2) Create a network consist of two subnetworks */ 386 /*-------------------------------------------------*/ 387 ierr = PetscLogStageRegister("Net Setup",&stage[1]);CHKERRQ(ierr); 388 ierr = PetscLogStagePush(stage[1]);CHKERRQ(ierr); 389 390 ierr = PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);CHKERRQ(ierr); 391 ierr = PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);CHKERRQ(ierr); 392 ierr = PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);CHKERRQ(ierr); 393 394 /* Create an empty network object */ 395 ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr); 396 397 /* Register the components in the network */ 398 ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);CHKERRQ(ierr); 399 ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);CHKERRQ(ierr); 400 ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);CHKERRQ(ierr); 401 ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);CHKERRQ(ierr); 402 403 ierr = DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);CHKERRQ(ierr); 404 ierr = DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);CHKERRQ(ierr); 405 #if 0 406 ierr = PetscPrintf(PETSC_COMM_WORLD,"power->compkey_branch %d\n",appctx_power->compkey_branch);CHKERRQ(ierr); 407 ierr = PetscPrintf(PETSC_COMM_WORLD,"power->compkey_bus %d\n",appctx_power->compkey_bus);CHKERRQ(ierr); 408 ierr = PetscPrintf(PETSC_COMM_WORLD,"power->compkey_gen %d\n",appctx_power->compkey_gen);CHKERRQ(ierr); 409 ierr = PetscPrintf(PETSC_COMM_WORLD,"power->compkey_load %d\n",appctx_power->compkey_load);CHKERRQ(ierr); 410 ierr = PetscPrintf(PETSC_COMM_WORLD,"water->compkey_edge %d\n",appctx_water->compkey_edge);CHKERRQ(ierr); 411 ierr = PetscPrintf(PETSC_COMM_WORLD,"water->compkey_vtx %d\n",appctx_water->compkey_vtx);CHKERRQ(ierr); 412 #endif 413 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Total local nvertices %D + %D = %D, nedges %D + %D = %D\n",rank,numVertices[0],numVertices[1],numVertices[0]+numVertices[1],numEdges[0],numEdges[1],numEdges[0]+numEdges[1]);CHKERRQ(ierr); 414 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 415 416 ierr = DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,Nsubnet);CHKERRQ(ierr); 417 ierr = DMNetworkAddSubnetwork(networkdm,"power",numEdges[0],edgelist_power,&power_netnum);CHKERRQ(ierr); 418 ierr = DMNetworkAddSubnetwork(networkdm,"water",numEdges[1],edgelist_water,&water_netnum);CHKERRQ(ierr); 419 420 /* vertex subnet[0].4 shares with vertex subnet[1].0 */ 421 power_svtx = 4; water_svtx = 0; 422 ierr = DMNetworkAddSharedVertices(networkdm,power_netnum,water_netnum,1,&power_svtx,&water_svtx);CHKERRQ(ierr); 423 424 /* Set up the network layout */ 425 ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr); 426 427 /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */ 428 /*-------------------------------------------------------*/ 429 genj = 0; loadj = 0; 430 ierr = DMNetworkGetSubnetwork(networkdm,power_netnum,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 431 432 for (i = 0; i < ne; i++) { 433 ierr = DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i],0);CHKERRQ(ierr); 434 } 435 436 for (i = 0; i < nv; i++) { 437 ierr = DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);CHKERRQ(ierr); 438 if (flg) continue; 439 440 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i],2);CHKERRQ(ierr); 441 if (pfdata->bus[i].ngen) { 442 for (j = 0; j < pfdata->bus[i].ngen; j++) { 443 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++],0);CHKERRQ(ierr); 444 } 445 } 446 if (pfdata->bus[i].nload) { 447 for (j=0; j < pfdata->bus[i].nload; j++) { 448 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++],0);CHKERRQ(ierr); 449 } 450 } 451 } 452 453 /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */ 454 /*-------------------------------------------------------*/ 455 ierr = DMNetworkGetSubnetwork(networkdm,water_netnum,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 456 for (i = 0; i < ne; i++) { 457 ierr = DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i],0);CHKERRQ(ierr); 458 } 459 460 for (i = 0; i < nv; i++) { 461 ierr = DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);CHKERRQ(ierr); 462 if (flg) continue; 463 464 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[i],1);CHKERRQ(ierr); 465 } 466 467 /* ADD VARIABLES AND COMPONENTS AT THE SHARED VERTEX: net[0].4 coupls with net[1].0 -- owning and all ghost ranks of the vertex do this */ 468 /*----------------------------------------------------------------------------------------------------------------------------*/ 469 ierr = DMNetworkGetSharedVertices(networkdm,&nv,&vtx);CHKERRQ(ierr); 470 for (i = 0; i < nv; i++) { 471 /* power */ 472 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[4],2);CHKERRQ(ierr); 473 /* bus[4] is a load, add its component */ 474 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[0],0);CHKERRQ(ierr); 475 476 /* water */ 477 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[0],1);CHKERRQ(ierr); 478 } 479 480 /* Set up DM for use */ 481 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 482 if (viewDM) { 483 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");CHKERRQ(ierr); 484 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 485 } 486 487 /* Free user objects */ 488 ierr = PetscFree(edgelist_power);CHKERRQ(ierr); 489 ierr = PetscFree(pfdata->bus);CHKERRQ(ierr); 490 ierr = PetscFree(pfdata->gen);CHKERRQ(ierr); 491 ierr = PetscFree(pfdata->branch);CHKERRQ(ierr); 492 ierr = PetscFree(pfdata->load);CHKERRQ(ierr); 493 ierr = PetscFree(pfdata);CHKERRQ(ierr); 494 495 ierr = PetscFree(edgelist_water);CHKERRQ(ierr); 496 ierr = PetscFree(waterdata->vertex);CHKERRQ(ierr); 497 ierr = PetscFree(waterdata->edge);CHKERRQ(ierr); 498 ierr = PetscFree(waterdata);CHKERRQ(ierr); 499 500 /* Re-distribute networkdm to multiple processes for better job balance */ 501 if (size >1 && distribute) { 502 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 503 if (viewDM) { 504 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");CHKERRQ(ierr); 505 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 506 } 507 } 508 509 /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */ 510 if (test) { 511 PetscInt v,gidx; 512 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 513 for (i=0; i<Nsubnet; i++) { 514 ierr = DMNetworkGetSubnetwork(networkdm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 515 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, subnet[%d] ne %d, nv %d\n",rank,i,ne,nv);CHKERRQ(ierr); 516 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 517 518 for (v=0; v<nv; v++) { 519 ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghost);CHKERRQ(ierr); 520 ierr = DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);CHKERRQ(ierr); 521 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] subnet[%d] v %d %d; ghost %d\n",rank,i,vtx[v],gidx,ghost);CHKERRQ(ierr); 522 } 523 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 524 } 525 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 526 527 ierr = DMNetworkGetSharedVertices(networkdm,&nv,&vtx);CHKERRQ(ierr); 528 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, num of shared vertices nsv = %d\n",rank,nv);CHKERRQ(ierr); 529 for (v=0; v<nv; v++) { 530 ierr = DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);CHKERRQ(ierr); 531 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] sv %d, gidx=%d\n",rank,vtx[v],gidx);CHKERRQ(ierr); 532 } 533 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 534 } 535 536 /* Create solution vector X */ 537 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 538 ierr = VecDuplicate(X,&F);CHKERRQ(ierr); 539 ierr = DMGetLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 540 PetscLogStagePop(); 541 542 /* (3) Setup Solvers */ 543 /*-------------------*/ 544 ierr = PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);CHKERRQ(ierr); 545 ierr = PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);CHKERRQ(ierr); 546 547 ierr = PetscLogStageRegister("SNES Setup",&stage[2]);CHKERRQ(ierr); 548 ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr); 549 550 ierr = SetInitialGuess(networkdm,X,&user);CHKERRQ(ierr); 551 552 /* Create coupled snes */ 553 /*-------------------- */ 554 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");CHKERRQ(ierr); 555 user.subsnes_id = Nsubnet; 556 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 557 ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr); 558 ierr = SNESSetOptionsPrefix(snes,"coupled_");CHKERRQ(ierr); 559 ierr = SNESSetFunction(snes,F,FormFunction,&user);CHKERRQ(ierr); 560 ierr = SNESMonitorSet(snes,UserMonitor,&user,NULL);CHKERRQ(ierr); 561 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 562 563 if (viewJ) { 564 /* View Jac structure */ 565 ierr = SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);CHKERRQ(ierr); 566 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 567 } 568 569 if (viewX) { 570 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");CHKERRQ(ierr); 571 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 572 } 573 574 if (viewJ) { 575 /* View assembled Jac */ 576 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 577 } 578 579 /* Create snes_power */ 580 /*-------------------*/ 581 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");CHKERRQ(ierr); 582 583 user.subsnes_id = 0; 584 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_power);CHKERRQ(ierr); 585 ierr = SNESSetDM(snes_power,networkdm);CHKERRQ(ierr); 586 ierr = SNESSetOptionsPrefix(snes_power,"power_");CHKERRQ(ierr); 587 ierr = SNESSetFunction(snes_power,F,FormFunction,&user);CHKERRQ(ierr); 588 ierr = SNESMonitorSet(snes_power,UserMonitor,&user,NULL);CHKERRQ(ierr); 589 590 /* Use user-provide Jacobian */ 591 ierr = DMCreateMatrix(networkdm,&Jac);CHKERRQ(ierr); 592 ierr = SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);CHKERRQ(ierr); 593 ierr = SNESSetFromOptions(snes_power);CHKERRQ(ierr); 594 595 if (viewX) { 596 ierr = PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");CHKERRQ(ierr); 597 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 598 } 599 if (viewJ) { 600 ierr = PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");CHKERRQ(ierr); 601 ierr = SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);CHKERRQ(ierr); 602 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 603 /* ierr = MatView(Jac,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 604 } 605 606 /* Create snes_water */ 607 /*-------------------*/ 608 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");CHKERRQ(ierr); 609 610 user.subsnes_id = 1; 611 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_water);CHKERRQ(ierr); 612 ierr = SNESSetDM(snes_water,networkdm);CHKERRQ(ierr); 613 ierr = SNESSetOptionsPrefix(snes_water,"water_");CHKERRQ(ierr); 614 ierr = SNESSetFunction(snes_water,F,FormFunction,&user);CHKERRQ(ierr); 615 ierr = SNESMonitorSet(snes_water,UserMonitor,&user,NULL);CHKERRQ(ierr); 616 ierr = SNESSetFromOptions(snes_water);CHKERRQ(ierr); 617 618 if (viewX) { 619 ierr = PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");CHKERRQ(ierr); 620 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 621 } 622 ierr = PetscLogStagePop();CHKERRQ(ierr); 623 624 /* (4) Solve */ 625 /*-----------*/ 626 ierr = PetscLogStageRegister("SNES Solve",&stage[3]);CHKERRQ(ierr); 627 ierr = PetscLogStagePush(stage[3]);CHKERRQ(ierr); 628 user.it = 0; 629 reason = SNES_DIVERGED_DTOL; 630 while (user.it < it_max && (PetscInt)reason<0) { 631 #if 0 632 user.subsnes_id = 0; 633 ierr = SNESSolve(snes_power,NULL,X);CHKERRQ(ierr); 634 635 user.subsnes_id = 1; 636 ierr = SNESSolve(snes_water,NULL,X);CHKERRQ(ierr); 637 #endif 638 user.subsnes_id = Nsubnet; 639 ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); 640 641 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 642 user.it++; 643 } 644 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);CHKERRQ(ierr); 645 if (viewX) { 646 ierr = PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");CHKERRQ(ierr); 647 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 648 } 649 ierr = PetscLogStagePop();CHKERRQ(ierr); 650 651 /* Free objects */ 652 /* -------------*/ 653 ierr = VecDestroy(&X);CHKERRQ(ierr); 654 ierr = VecDestroy(&F);CHKERRQ(ierr); 655 ierr = DMRestoreLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 656 657 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 658 ierr = MatDestroy(&Jac);CHKERRQ(ierr); 659 ierr = SNESDestroy(&snes_power);CHKERRQ(ierr); 660 ierr = SNESDestroy(&snes_water);CHKERRQ(ierr); 661 662 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 663 ierr = PetscFinalize(); 664 return ierr; 665 } 666 667 /*TEST 668 669 build: 670 requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED) 671 depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c 672 673 test: 674 args: -coupled_snes_converged_reason -options_left no -viewDM 675 localrunfiles: ex1options power/case9.m water/sample1.inp 676 output_file: output/ex1.out 677 678 test: 679 suffix: 2 680 nsize: 3 681 args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis 682 localrunfiles: ex1options power/case9.m water/sample1.inp 683 output_file: output/ex1_2.out 684 requires: parmetis 685 686 # test: 687 # suffix: 3 688 # nsize: 3 689 # args: -coupled_snes_converged_reason -options_left no -distribute false 690 # localrunfiles: ex1options power/case9.m water/sample1.inp 691 # output_file: output/ex1_2.out 692 693 test: 694 suffix: 4 695 nsize: 4 696 args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -viewDM 697 localrunfiles: ex1options power/case9.m water/sample1.inp 698 output_file: output/ex1_4.out 699 700 TEST*/ 701