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) { 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) 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) { 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",numVertices[0],numEdges[0],edgelist_power,&power_netnum);CHKERRQ(ierr); 418 ierr = DMNetworkAddSubnetwork(networkdm,"water",numVertices[1],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 -- only the owner of the vertex does this */ 468 /*----------------------------------------------------------------------------------------------------------------------------*/ 469 ierr = DMNetworkGetSharedVertices(networkdm,&nv,&vtx);CHKERRQ(ierr); 470 for (i = 0; i < nv; i++) { 471 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);CHKERRQ(ierr); 472 /* printf("[%d] coupling info: nv %d; sv[0] %d; ghost %d\n",rank,nv,vtx[0],ghost);CHKERRQ(ierr); */ 473 if (ghost) continue; 474 475 /* power */ 476 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[4],2);CHKERRQ(ierr); 477 /* bus[4] is a load, add its component */ 478 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[0],0);CHKERRQ(ierr); 479 480 /* water */ 481 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[0],1);CHKERRQ(ierr); 482 } 483 484 /* Set up DM for use */ 485 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 486 if (viewDM) { 487 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");CHKERRQ(ierr); 488 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 489 } 490 491 /* Free user objects */ 492 ierr = PetscFree(edgelist_power);CHKERRQ(ierr); 493 ierr = PetscFree(pfdata->bus);CHKERRQ(ierr); 494 ierr = PetscFree(pfdata->gen);CHKERRQ(ierr); 495 ierr = PetscFree(pfdata->branch);CHKERRQ(ierr); 496 ierr = PetscFree(pfdata->load);CHKERRQ(ierr); 497 ierr = PetscFree(pfdata);CHKERRQ(ierr); 498 499 ierr = PetscFree(edgelist_water);CHKERRQ(ierr); 500 ierr = PetscFree(waterdata->vertex);CHKERRQ(ierr); 501 ierr = PetscFree(waterdata->edge);CHKERRQ(ierr); 502 ierr = PetscFree(waterdata);CHKERRQ(ierr); 503 504 /* Re-distribute networkdm to multiple processes for better job balance */ 505 if (size >1 && distribute) { 506 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 507 if (viewDM) { 508 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");CHKERRQ(ierr); 509 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 510 } 511 } 512 513 /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */ 514 if (test) { 515 PetscInt v,gidx; 516 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 517 for (i=0; i<Nsubnet; i++) { 518 ierr = DMNetworkGetSubnetwork(networkdm,i,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 519 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, subnet[%d] ne %d, nv %d\n",rank,i,ne,nv);CHKERRQ(ierr); 520 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 521 522 for (v=0; v<nv; v++) { 523 ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghost);CHKERRQ(ierr); 524 ierr = DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);CHKERRQ(ierr); 525 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] subnet[%d] v %d %d; ghost %d\n",rank,i,vtx[v],gidx,ghost);CHKERRQ(ierr); 526 } 527 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 528 } 529 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 530 531 ierr = DMNetworkGetSharedVertices(networkdm,&nv,&vtx);CHKERRQ(ierr); 532 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, num of shared vertices nsv = %d\n",rank,nv);CHKERRQ(ierr); 533 for (v=0; v<nv; v++) { 534 ierr = DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);CHKERRQ(ierr); 535 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] sv %d, gidx=%d\n",rank,vtx[v],gidx);CHKERRQ(ierr); 536 } 537 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 538 } 539 540 /* Create solution vector X */ 541 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 542 ierr = VecDuplicate(X,&F);CHKERRQ(ierr); 543 ierr = DMGetLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 544 PetscLogStagePop(); 545 546 /* (3) Setup Solvers */ 547 /*-------------------*/ 548 ierr = PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);CHKERRQ(ierr); 549 ierr = PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);CHKERRQ(ierr); 550 551 ierr = PetscLogStageRegister("SNES Setup",&stage[2]);CHKERRQ(ierr); 552 ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr); 553 554 ierr = SetInitialGuess(networkdm,X,&user);CHKERRQ(ierr); 555 556 /* Create coupled snes */ 557 /*-------------------- */ 558 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");CHKERRQ(ierr); 559 user.subsnes_id = Nsubnet; 560 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 561 ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr); 562 ierr = SNESSetOptionsPrefix(snes,"coupled_");CHKERRQ(ierr); 563 ierr = SNESSetFunction(snes,F,FormFunction,&user);CHKERRQ(ierr); 564 ierr = SNESMonitorSet(snes,UserMonitor,&user,NULL);CHKERRQ(ierr); 565 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 566 567 if (viewJ) { 568 /* View Jac structure */ 569 ierr = SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);CHKERRQ(ierr); 570 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 571 } 572 573 if (viewX) { 574 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");CHKERRQ(ierr); 575 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 576 } 577 578 if (viewJ) { 579 /* View assembled Jac */ 580 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 581 } 582 583 /* Create snes_power */ 584 /*-------------------*/ 585 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");CHKERRQ(ierr); 586 587 user.subsnes_id = 0; 588 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_power);CHKERRQ(ierr); 589 ierr = SNESSetDM(snes_power,networkdm);CHKERRQ(ierr); 590 ierr = SNESSetOptionsPrefix(snes_power,"power_");CHKERRQ(ierr); 591 ierr = SNESSetFunction(snes_power,F,FormFunction,&user);CHKERRQ(ierr); 592 ierr = SNESMonitorSet(snes_power,UserMonitor,&user,NULL);CHKERRQ(ierr); 593 594 /* Use user-provide Jacobian */ 595 ierr = DMCreateMatrix(networkdm,&Jac);CHKERRQ(ierr); 596 ierr = SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);CHKERRQ(ierr); 597 ierr = SNESSetFromOptions(snes_power);CHKERRQ(ierr); 598 599 if (viewX) { 600 ierr = PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");CHKERRQ(ierr); 601 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 602 } 603 if (viewJ) { 604 ierr = PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");CHKERRQ(ierr); 605 ierr = SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);CHKERRQ(ierr); 606 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 607 /* ierr = MatView(Jac,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 608 } 609 610 /* Create snes_water */ 611 /*-------------------*/ 612 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");CHKERRQ(ierr); 613 614 user.subsnes_id = 1; 615 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_water);CHKERRQ(ierr); 616 ierr = SNESSetDM(snes_water,networkdm);CHKERRQ(ierr); 617 ierr = SNESSetOptionsPrefix(snes_water,"water_");CHKERRQ(ierr); 618 ierr = SNESSetFunction(snes_water,F,FormFunction,&user);CHKERRQ(ierr); 619 ierr = SNESMonitorSet(snes_water,UserMonitor,&user,NULL);CHKERRQ(ierr); 620 ierr = SNESSetFromOptions(snes_water);CHKERRQ(ierr); 621 622 if (viewX) { 623 ierr = PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");CHKERRQ(ierr); 624 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 625 } 626 ierr = PetscLogStagePop();CHKERRQ(ierr); 627 628 /* (4) Solve */ 629 /*-----------*/ 630 ierr = PetscLogStageRegister("SNES Solve",&stage[3]);CHKERRQ(ierr); 631 ierr = PetscLogStagePush(stage[3]);CHKERRQ(ierr); 632 user.it = 0; 633 reason = SNES_DIVERGED_DTOL; 634 while (user.it < it_max && (PetscInt)reason<0) { 635 #if 0 636 user.subsnes_id = 0; 637 ierr = SNESSolve(snes_power,NULL,X);CHKERRQ(ierr); 638 639 user.subsnes_id = 1; 640 ierr = SNESSolve(snes_water,NULL,X);CHKERRQ(ierr); 641 #endif 642 user.subsnes_id = Nsubnet; 643 ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); 644 645 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 646 user.it++; 647 } 648 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);CHKERRQ(ierr); 649 if (viewX) { 650 ierr = PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");CHKERRQ(ierr); 651 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 652 } 653 ierr = PetscLogStagePop();CHKERRQ(ierr); 654 655 /* Free objects */ 656 /* -------------*/ 657 ierr = VecDestroy(&X);CHKERRQ(ierr); 658 ierr = VecDestroy(&F);CHKERRQ(ierr); 659 ierr = DMRestoreLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 660 661 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 662 ierr = MatDestroy(&Jac);CHKERRQ(ierr); 663 ierr = SNESDestroy(&snes_power);CHKERRQ(ierr); 664 ierr = SNESDestroy(&snes_water);CHKERRQ(ierr); 665 666 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 667 ierr = PetscFinalize(); 668 return ierr; 669 } 670 671 /*TEST 672 673 build: 674 requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED) 675 depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c 676 677 test: 678 args: -coupled_snes_converged_reason -options_left no -viewDM 679 localrunfiles: ex1options power/case9.m water/sample1.inp 680 output_file: output/ex1.out 681 682 test: 683 suffix: 2 684 nsize: 3 685 args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis 686 localrunfiles: ex1options power/case9.m water/sample1.inp 687 output_file: output/ex1_2.out 688 requires: parmetis 689 690 # test: 691 # suffix: 3 692 # nsize: 3 693 # args: -coupled_snes_converged_reason -options_left no -distribute false 694 # localrunfiles: ex1options power/case9.m water/sample1.inp 695 # output_file: output/ex1_2.out 696 697 test: 698 suffix: 4 699 nsize: 4 700 args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -viewDM 701 localrunfiles: ex1options power/case9.m water/sample1.inp 702 output_file: output/ex1_4.out 703 704 TEST*/ 705