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);CHKERRQ(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,fnorm);CHKERRQ(ierr); 41 } else { 42 ierr = PetscPrintf(PETSC_COMM_SELF," subsnes_id %d, fnorm %g\n",user->subsnes_id,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);CHKERRQ(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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetVariableGlobalOffset(networkdm,vtx[i],&offset);CHKERRQ(ierr); 88 ierr = DMNetworkGetNumVariables(networkdm,vtx[i],&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 = DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);CHKERRQ(ierr); 122 ierr = DMNetworkGetNumVariables(networkdm,vtx[i],&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; 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 146 PetscFunctionBegin; 147 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 148 ierr = PetscObjectGetComm((PetscObject)networkdm,&comm);CHKERRQ(ierr); 149 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 150 151 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 152 ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr); 153 ierr = VecSet(F,0.0);CHKERRQ(ierr); 154 ierr = VecSet(localF,0.0);CHKERRQ(ierr); 155 156 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 157 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 158 159 /* Form Function for power subnetwork */ 160 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 161 if (user->subsnes_id == 1) { /* snes_water only */ 162 ierr = FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);CHKERRQ(ierr); 163 } else { 164 ierr = FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,&appctx_power);CHKERRQ(ierr); 165 } 166 167 /* Form Function for water subnetwork */ 168 ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 169 if (user->subsnes_id == 0) { /* snes_power only */ 170 ierr = FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);CHKERRQ(ierr); 171 } else { 172 ierr = FormFunction_Water(networkdm,localX,localF,nv,ne,vtx,edges,NULL);CHKERRQ(ierr); 173 } 174 175 #if 0 176 /* Form Function for the coupling subnetwork -- experimental, not done yet */ 177 ierr = DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);CHKERRQ(ierr); 178 if (ne) { 179 const PetscInt *cone,*connedges,*econe; 180 PetscInt key,vid[2],nconnedges,k,e,keye; 181 void* component; 182 AppCtx_Water appctx_water = (*user).appctx_water; 183 184 ierr = DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);CHKERRQ(ierr); 185 186 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);CHKERRQ(ierr); 187 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);CHKERRQ(ierr); 188 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Formfunction, coupling subnetwork: ne %d; connected vertices %d %d\n",rank,ne,vid[0],vid[1]);CHKERRQ(ierr); 189 190 /* Get coupling powernet load vertex */ 191 ierr = DMNetworkGetComponent(networkdm,cone[0],1,&key,&component);CHKERRQ(ierr); 192 if (key != appctx_power.compkey_load) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power load vertex"); 193 194 /* Get coupling water vertex and pump edge */ 195 ierr = DMNetworkGetComponent(networkdm,cone[1],0,&key,&component);CHKERRQ(ierr); 196 if (key != appctx_water.compkey_vtx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex"); 197 198 /* get its supporting edges */ 199 ierr = DMNetworkGetSupportingEdges(networkdm,cone[1],&nconnedges,&connedges);CHKERRQ(ierr); 200 201 for (k=0; k<nconnedges; k++) { 202 e = connedges[k]; 203 ierr = DMNetworkGetComponent(networkdm,e,0,&keye,&component);CHKERRQ(ierr); 204 205 if (keye == appctx_water.compkey_edge) { /* water_compkey_edge */ 206 EDGE_Water edge=(EDGE_Water)component; 207 if (edge->type == EDGE_TYPE_PUMP) { 208 /* compute flow_pump */ 209 PetscInt offsetnode1,offsetnode2,key_0,key_1; 210 const PetscScalar *xarr; 211 PetscScalar *farr; 212 VERTEX_Water vertexnode1,vertexnode2; 213 214 ierr = DMNetworkGetConnectedVertices(networkdm,e,&econe);CHKERRQ(ierr); 215 ierr = DMNetworkGetGlobalVertexIndex(networkdm,econe[0],&vid[0]);CHKERRQ(ierr); 216 ierr = DMNetworkGetGlobalVertexIndex(networkdm,econe[1],&vid[1]);CHKERRQ(ierr); 217 218 ierr = VecGetArray(localF,&farr);CHKERRQ(ierr); 219 ierr = DMNetworkGetVariableOffset(networkdm,econe[0],&offsetnode1);CHKERRQ(ierr); 220 ierr = DMNetworkGetVariableOffset(networkdm,econe[1],&offsetnode2);CHKERRQ(ierr); 221 222 ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); 223 #if 0 224 PetscScalar hf,ht; 225 Pump *pump; 226 pump = &edge->pump; 227 hf = xarr[offsetnode1]; 228 ht = xarr[offsetnode2]; 229 230 PetscScalar flow = Flow_Pump(pump,hf,ht); 231 PetscScalar Hp = 0.1; /* load->pl */ 232 PetscScalar flow_couple = 8.81*Hp*1.e6/(ht-hf); /* pump->h0; */ 233 /* ierr = PetscPrintf(PETSC_COMM_SELF,"pump %d: connected vtx %d %d; flow_pump %g flow_couple %g; offset %d %d\n",e,vid[0],vid[1],flow,flow_couple,offsetnode1,offsetnode2);CHKERRQ(ierr); */ 234 #endif 235 /* Get the components at the two vertices */ 236 ierr = DMNetworkGetComponent(networkdm,econe[0],0,&key_0,(void**)&vertexnode1);CHKERRQ(ierr); 237 if (key_0 != appctx_water.compkey_vtx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex"); 238 ierr = DMNetworkGetComponent(networkdm,econe[1],0,&key_1,(void**)&vertexnode2);CHKERRQ(ierr); 239 if (key_1 != appctx_water.compkey_vtx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex"); 240 #if 0 241 /* subtract flow_pump computed in FormFunction_Water() and add flow_couple to connected nodes */ 242 if (vertexnode1->type == VERTEX_TYPE_JUNCTION) { 243 farr[offsetnode1] += flow; 244 farr[offsetnode1] -= flow_couple; 245 } 246 if (vertexnode2->type == VERTEX_TYPE_JUNCTION) { 247 farr[offsetnode2] -= flow; 248 farr[offsetnode2] += flow_couple; 249 } 250 #endif 251 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 252 ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); 253 } 254 } 255 } 256 } 257 #endif 258 259 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 260 261 ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 262 ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 263 ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr); 264 /* ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 265 PetscFunctionReturn(0); 266 } 267 268 PetscErrorCode SetInitialGuess(DM networkdm,Vec X,void* appctx) 269 { 270 PetscErrorCode ierr; 271 PetscInt nv,ne; 272 const PetscInt *vtx,*edges; 273 UserCtx *user = (UserCtx*)appctx; 274 Vec localX = user->localXold; 275 UserCtx_Power appctx_power = (*user).appctx_power; 276 277 PetscFunctionBegin; 278 ierr = VecSet(X,0.0);CHKERRQ(ierr); 279 ierr = VecSet(localX,0.0);CHKERRQ(ierr); 280 281 /* Set initial guess for power subnetwork */ 282 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 283 ierr = SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,&appctx_power);CHKERRQ(ierr); 284 285 /* Set initial guess for water subnetwork */ 286 ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 287 ierr = SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);CHKERRQ(ierr); 288 289 #if 0 290 /* Set initial guess for the coupling subnet */ 291 ierr = DMNetworkGetSubnetworkInfo(networkdm,2,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 292 if (ne) { 293 const PetscInt *cone; 294 ierr = DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);CHKERRQ(ierr); 295 } 296 #endif 297 298 ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 299 ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 int main(int argc,char **argv) 304 { 305 PetscErrorCode ierr; 306 DM networkdm; 307 PetscLogStage stage[4]; 308 PetscMPIInt rank,size; 309 PetscInt nsubnet=2,nsubnetCouple=0,numVertices[2],numEdges[2],numEdgesCouple[1]; 310 PetscInt i,j,nv,ne; 311 PetscInt *edgelist[2]; 312 const PetscInt *vtx,*edges; 313 Vec X,F; 314 SNES snes,snes_power,snes_water; 315 Mat Jac; 316 PetscBool viewJ=PETSC_FALSE,viewX=PETSC_FALSE,viewDM=PETSC_FALSE,test=PETSC_FALSE,distribute=PETSC_TRUE; 317 UserCtx user; 318 PetscInt it_max=10; 319 SNESConvergedReason reason; 320 321 /* Power subnetwork */ 322 UserCtx_Power *appctx_power = &user.appctx_power; 323 char pfdata_file[PETSC_MAX_PATH_LEN]="power/case9.m"; 324 PFDATA *pfdata=NULL; 325 PetscInt genj,loadj; 326 PetscInt *edgelist_power=NULL; 327 PetscScalar Sbase=0.0; 328 329 /* Water subnetwork */ 330 AppCtx_Water *appctx_water = &user.appctx_water; 331 WATERDATA *waterdata=NULL; 332 char waterdata_file[PETSC_MAX_PATH_LEN]="water/sample1.inp"; 333 PetscInt *edgelist_water=NULL; 334 335 /* Coupling subnetwork */ 336 PetscInt *edgelist_couple=NULL; 337 338 ierr = PetscInitialize(&argc,&argv,"ex1options",help);if (ierr) return ierr; 339 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 340 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 341 342 /* (1) Read Data - Only rank 0 reads the data */ 343 /*--------------------------------------------*/ 344 ierr = PetscLogStageRegister("Read Data",&stage[0]);CHKERRQ(ierr); 345 PetscLogStagePush(stage[0]); 346 347 for (i=0; i<nsubnet; i++) { 348 numVertices[i] = 0; 349 numEdges[i] = 0; 350 } 351 numEdgesCouple[0] = 0; 352 353 /* proc[0] READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */ 354 if (rank == 0) { 355 ierr = PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);CHKERRQ(ierr); 356 ierr = PetscNew(&pfdata);CHKERRQ(ierr); 357 ierr = PFReadMatPowerData(pfdata,pfdata_file);CHKERRQ(ierr); 358 Sbase = pfdata->sbase; 359 360 numEdges[0] = pfdata->nbranch; 361 numVertices[0] = pfdata->nbus; 362 363 ierr = PetscMalloc1(2*numEdges[0],&edgelist_power);CHKERRQ(ierr); 364 ierr = GetListofEdges_Power(pfdata,edgelist_power);CHKERRQ(ierr); 365 } 366 /* Broadcast power Sbase to all processors */ 367 ierr = MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr); 368 appctx_power->Sbase = Sbase; 369 appctx_power->jac_error = PETSC_FALSE; 370 /* If external option activated. Introduce error in jacobian */ 371 ierr = PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);CHKERRQ(ierr); 372 373 /* proc[1] GET DATA FOR THE SECOND SUBNETWORK: Water */ 374 if (size == 1 || (size > 1 && rank == 1)) { 375 ierr = PetscNew(&waterdata);CHKERRQ(ierr); 376 ierr = PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,sizeof(waterdata_file),NULL);CHKERRQ(ierr); 377 ierr = WaterReadData(waterdata,waterdata_file);CHKERRQ(ierr); 378 379 ierr = PetscCalloc1(2*waterdata->nedge,&edgelist_water);CHKERRQ(ierr); 380 ierr = GetListofEdges_Water(waterdata,edgelist_water);CHKERRQ(ierr); 381 numEdges[1] = waterdata->nedge; 382 numVertices[1] = waterdata->nvertex; 383 } 384 385 /* Get data for the coupling subnetwork */ 386 if (size == 1) { /* TODO: for size > 1, parallel processing coupling is buggy */ 387 nsubnetCouple = 1; 388 numEdgesCouple[0] = 1; 389 390 ierr = PetscMalloc1(4*numEdgesCouple[0],&edgelist_couple);CHKERRQ(ierr); 391 edgelist_couple[0] = 0; edgelist_couple[1] = 4; /* from node: net[0] vertex[4] */ 392 edgelist_couple[2] = 1; edgelist_couple[3] = 0; /* to node: net[1] vertex[0] */ 393 } 394 PetscLogStagePop(); 395 396 /* (2) Create network */ 397 /*--------------------*/ 398 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr); 399 ierr = PetscLogStageRegister("Net Setup",&stage[1]);CHKERRQ(ierr); 400 PetscLogStagePush(stage[1]); 401 402 ierr = PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);CHKERRQ(ierr); 403 ierr = PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);CHKERRQ(ierr); 404 ierr = PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);CHKERRQ(ierr); 405 406 /* Create an empty network object */ 407 ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr); 408 409 /* Register the components in the network */ 410 ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);CHKERRQ(ierr); 411 ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);CHKERRQ(ierr); 412 ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);CHKERRQ(ierr); 413 ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);CHKERRQ(ierr); 414 415 ierr = DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);CHKERRQ(ierr); 416 ierr = DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);CHKERRQ(ierr); 417 418 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Total local nvertices %D + %D = %D, nedges %D + %D + %D = %D\n",rank,numVertices[0],numVertices[1],numVertices[0]+numVertices[1],numEdges[0],numEdges[1],numEdgesCouple[0],numEdges[0]+numEdges[1]+numEdgesCouple[0]);CHKERRQ(ierr); 419 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 420 421 ierr = DMNetworkSetSizes(networkdm,nsubnet,numVertices,numEdges,nsubnetCouple,numEdgesCouple);CHKERRQ(ierr); 422 423 /* Add edge connectivity */ 424 edgelist[0] = edgelist_power; 425 edgelist[1] = edgelist_water; 426 ierr = DMNetworkSetEdgeList(networkdm,edgelist,&edgelist_couple);CHKERRQ(ierr); 427 428 /* Set up the network layout */ 429 ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr); 430 431 /* Add network components - only process[0] has any data to add */ 432 /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */ 433 genj = 0; loadj = 0; 434 if (rank == 0) { 435 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 436 /* ierr = PetscPrintf(PETSC_COMM_SELF,"Power network: nv %D, ne %D\n",nv,ne);CHKERRQ(ierr); */ 437 for (i = 0; i < ne; i++) { 438 ierr = DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i]);CHKERRQ(ierr); 439 } 440 441 for (i = 0; i < nv; i++) { 442 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i]);CHKERRQ(ierr); 443 if (pfdata->bus[i].ngen) { 444 for (j = 0; j < pfdata->bus[i].ngen; j++) { 445 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++]);CHKERRQ(ierr); 446 } 447 } 448 if (pfdata->bus[i].nload) { 449 for (j=0; j < pfdata->bus[i].nload; j++) { 450 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++]);CHKERRQ(ierr); 451 } 452 } 453 /* Add number of variables */ 454 ierr = DMNetworkAddNumVariables(networkdm,vtx[i],2);CHKERRQ(ierr); 455 } 456 } 457 458 /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */ 459 if (size == 1 || (size > 1 && rank == 1)) { 460 ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 461 /* ierr = PetscPrintf(PETSC_COMM_SELF,"Water network: nv %D, ne %D\n",nv,ne);CHKERRQ(ierr); */ 462 for (i = 0; i < ne; i++) { 463 ierr = DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i]);CHKERRQ(ierr); 464 } 465 466 for (i = 0; i < nv; i++) { 467 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[i]);CHKERRQ(ierr); 468 /* Add number of variables */ 469 ierr = DMNetworkAddNumVariables(networkdm,vtx[i],1);CHKERRQ(ierr); 470 } 471 } 472 473 /* Set up DM for use */ 474 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 475 if (viewDM) { 476 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");CHKERRQ(ierr); 477 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 478 } 479 480 ierr = DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);CHKERRQ(ierr); 481 if (ne && viewDM) { 482 const PetscInt *cone; 483 PetscInt vid[2]; 484 ierr = DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);CHKERRQ(ierr); 485 486 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);CHKERRQ(ierr); 487 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);CHKERRQ(ierr); 488 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Coupling subnetwork: ne %d; connected vertices %d %d\n",rank,ne,vid[0],vid[1]);CHKERRQ(ierr); 489 } 490 491 /* Free user objects */ 492 if (!rank) { 493 ierr = PetscFree(edgelist_power);CHKERRQ(ierr); 494 ierr = PetscFree(pfdata->bus);CHKERRQ(ierr); 495 ierr = PetscFree(pfdata->gen);CHKERRQ(ierr); 496 ierr = PetscFree(pfdata->branch);CHKERRQ(ierr); 497 ierr = PetscFree(pfdata->load);CHKERRQ(ierr); 498 ierr = PetscFree(pfdata);CHKERRQ(ierr); 499 } 500 if (size == 1 || (size > 1 && rank == 1)) { 501 ierr = PetscFree(edgelist_water);CHKERRQ(ierr); 502 ierr = PetscFree(waterdata->vertex);CHKERRQ(ierr); 503 ierr = PetscFree(waterdata->edge);CHKERRQ(ierr); 504 ierr = PetscFree(waterdata);CHKERRQ(ierr); 505 } 506 if (size == 1) { 507 ierr = PetscFree(edgelist_couple);CHKERRQ(ierr); 508 } 509 510 /* Re-distribute networkdm to multiple processes for better job balance */ 511 if (distribute) { 512 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 513 if (viewDM) { 514 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");CHKERRQ(ierr); 515 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 516 } 517 } 518 519 /* Test DMNetworkGetSubnetworkCoupleInfo() */ 520 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr); 521 if (test) { 522 for (i=0; i<nsubnetCouple; i++) { 523 ierr = DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);CHKERRQ(ierr); 524 if (ne && viewDM) { 525 const PetscInt *cone; 526 PetscInt vid[2]; 527 ierr = DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);CHKERRQ(ierr); 528 529 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);CHKERRQ(ierr); 530 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);CHKERRQ(ierr); 531 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] After DMNetworkDistribute(), subnetworkCouple[%d]: ne %d; connected vertices %d %d\n",rank,i,ne,vid[0],vid[1]);CHKERRQ(ierr); 532 } 533 } 534 } 535 536 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 537 ierr = VecDuplicate(X,&F);CHKERRQ(ierr); 538 ierr = DMGetLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 539 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 PetscLogStagePush(stage[2]); 549 550 ierr = SetInitialGuess(networkdm,X,&user);CHKERRQ(ierr); 551 /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 552 553 /* Create coupled snes */ 554 /*-------------------- */ 555 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");CHKERRQ(ierr); 556 user.subsnes_id = nsubnet; 557 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 558 ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr); 559 ierr = SNESSetOptionsPrefix(snes,"coupled_");CHKERRQ(ierr); 560 ierr = SNESSetFunction(snes,F,FormFunction,&user);CHKERRQ(ierr); 561 ierr = SNESMonitorSet(snes,UserMonitor,&user,NULL);CHKERRQ(ierr); 562 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 563 564 if (viewJ) { 565 /* View Jac structure */ 566 ierr = SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);CHKERRQ(ierr); 567 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 568 } 569 570 if (viewX) { 571 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");CHKERRQ(ierr); 572 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 573 } 574 575 if (viewJ) { 576 /* View assembled Jac */ 577 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 578 } 579 580 /* Create snes_power */ 581 /*-------------------*/ 582 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");CHKERRQ(ierr); 583 ierr = SetInitialGuess(networkdm,X,&user);CHKERRQ(ierr); 584 585 user.subsnes_id = 0; 586 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_power);CHKERRQ(ierr); 587 ierr = SNESSetDM(snes_power,networkdm);CHKERRQ(ierr); 588 ierr = SNESSetOptionsPrefix(snes_power,"power_");CHKERRQ(ierr); 589 ierr = SNESSetFunction(snes_power,F,FormFunction,&user);CHKERRQ(ierr); 590 ierr = SNESMonitorSet(snes_power,UserMonitor,&user,NULL);CHKERRQ(ierr); 591 592 /* Use user-provide Jacobian */ 593 ierr = DMCreateMatrix(networkdm,&Jac);CHKERRQ(ierr); 594 ierr = SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);CHKERRQ(ierr); 595 ierr = SNESSetFromOptions(snes_power);CHKERRQ(ierr); 596 597 if (viewX) { 598 ierr = PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");CHKERRQ(ierr); 599 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 600 } 601 if (viewJ) { 602 ierr = PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");CHKERRQ(ierr); 603 ierr = SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);CHKERRQ(ierr); 604 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 605 /* ierr = MatView(Jac,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 606 } 607 608 /* Create snes_water */ 609 /*-------------------*/ 610 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");CHKERRQ(ierr); 611 ierr = SetInitialGuess(networkdm,X,&user);CHKERRQ(ierr); 612 613 user.subsnes_id = 1; 614 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_water);CHKERRQ(ierr); 615 ierr = SNESSetDM(snes_water,networkdm);CHKERRQ(ierr); 616 ierr = SNESSetOptionsPrefix(snes_water,"water_");CHKERRQ(ierr); 617 ierr = SNESSetFunction(snes_water,F,FormFunction,&user);CHKERRQ(ierr); 618 ierr = SNESMonitorSet(snes_water,UserMonitor,&user,NULL);CHKERRQ(ierr); 619 ierr = SNESSetFromOptions(snes_water);CHKERRQ(ierr); 620 621 if (viewX) { 622 ierr = PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");CHKERRQ(ierr); 623 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 624 } 625 PetscLogStagePop(); 626 627 /* (4) Solve */ 628 /*-----------*/ 629 ierr = PetscLogStageRegister("SNES Solve",&stage[3]);CHKERRQ(ierr); 630 PetscLogStagePush(stage[3]); 631 user.it = 0; 632 reason = SNES_DIVERGED_DTOL; 633 while (user.it < it_max && (PetscInt)reason<0) { 634 #if 0 635 user.subsnes_id = 0; 636 ierr = SNESSolve(snes_power,NULL,X);CHKERRQ(ierr); 637 638 user.subsnes_id = 1; 639 ierr = SNESSolve(snes_water,NULL,X);CHKERRQ(ierr); 640 #endif 641 user.subsnes_id = nsubnet; 642 ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); 643 644 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 645 user.it++; 646 } 647 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);CHKERRQ(ierr); 648 if (viewX) { 649 ierr = PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");CHKERRQ(ierr); 650 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 651 } 652 PetscLogStagePop(); 653 654 /* Free objects */ 655 /* -------------*/ 656 ierr = VecDestroy(&X);CHKERRQ(ierr); 657 ierr = VecDestroy(&F);CHKERRQ(ierr); 658 ierr = DMRestoreLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 659 660 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 661 ierr = MatDestroy(&Jac);CHKERRQ(ierr); 662 ierr = SNESDestroy(&snes_power);CHKERRQ(ierr); 663 ierr = SNESDestroy(&snes_water);CHKERRQ(ierr); 664 665 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 666 ierr = PetscFinalize(); 667 return ierr; 668 } 669 670 /*TEST 671 672 build: 673 requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED) 674 depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c 675 676 test: 677 args: -coupled_snes_converged_reason -options_left no 678 localrunfiles: ex1options power/case9.m water/sample1.inp 679 output_file: output/ex1.out 680 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 681 682 test: 683 suffix: 2 684 nsize: 3 685 args: -coupled_snes_converged_reason -options_left no 686 localrunfiles: ex1options power/case9.m water/sample1.inp 687 output_file: output/ex1_2.out 688 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 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 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 697 698 TEST*/ 699