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,(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);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 #if defined(PETSC_USE_LOG) 308 PetscLogStage stage[4]; 309 #endif 310 PetscMPIInt rank,size; 311 PetscInt nsubnet = 2,nsubnetCouple = 0,numVertices[2],numEdges[2],numEdgesCouple[1]; 312 PetscInt i,j,nv,ne; 313 PetscInt *edgelist[2]; 314 const PetscInt *vtx,*edges; 315 Vec X,F; 316 SNES snes,snes_power,snes_water; 317 Mat Jac; 318 PetscBool viewJ = PETSC_FALSE,viewX = PETSC_FALSE,viewDM = PETSC_FALSE,test = PETSC_FALSE,distribute = PETSC_TRUE; 319 UserCtx user; 320 PetscInt it_max = 10; 321 SNESConvergedReason reason; 322 323 /* Power subnetwork */ 324 UserCtx_Power *appctx_power = &user.appctx_power; 325 char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m"; 326 PFDATA *pfdata = NULL; 327 PetscInt genj,loadj; 328 PetscInt *edgelist_power = NULL; 329 PetscScalar Sbase = 0.0; 330 331 /* Water subnetwork */ 332 AppCtx_Water *appctx_water = &user.appctx_water; 333 WATERDATA *waterdata = NULL; 334 char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp"; 335 PetscInt *edgelist_water = NULL; 336 337 /* Coupling subnetwork */ 338 PetscInt *edgelist_couple = NULL; 339 340 ierr = PetscInitialize(&argc,&argv,"ex1options",help);if (ierr) return ierr; 341 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 342 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 343 344 /* (1) Read Data - Only rank 0 reads the data */ 345 /*--------------------------------------------*/ 346 ierr = PetscLogStageRegister("Read Data",&stage[0]);CHKERRQ(ierr); 347 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 348 349 for (i=0; i<nsubnet; i++) { 350 numVertices[i] = 0; 351 numEdges[i] = 0; 352 } 353 numEdgesCouple[0] = 0; 354 355 /* proc[0] READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */ 356 if (rank == 0) { 357 ierr = PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);CHKERRQ(ierr); 358 ierr = PetscNew(&pfdata);CHKERRQ(ierr); 359 ierr = PFReadMatPowerData(pfdata,pfdata_file);CHKERRQ(ierr); 360 Sbase = pfdata->sbase; 361 362 numEdges[0] = pfdata->nbranch; 363 numVertices[0] = pfdata->nbus; 364 365 ierr = PetscMalloc1(2*numEdges[0],&edgelist_power);CHKERRQ(ierr); 366 ierr = GetListofEdges_Power(pfdata,edgelist_power);CHKERRQ(ierr); 367 } 368 /* Broadcast power Sbase to all processors */ 369 ierr = MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr); 370 appctx_power->Sbase = Sbase; 371 appctx_power->jac_error = PETSC_FALSE; 372 /* If external option activated. Introduce error in jacobian */ 373 ierr = PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);CHKERRQ(ierr); 374 375 /* proc[1] GET DATA FOR THE SECOND SUBNETWORK: Water */ 376 if (size == 1 || (size > 1 && rank == 1)) { 377 ierr = PetscNew(&waterdata);CHKERRQ(ierr); 378 ierr = PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,sizeof(waterdata_file),NULL);CHKERRQ(ierr); 379 ierr = WaterReadData(waterdata,waterdata_file);CHKERRQ(ierr); 380 381 ierr = PetscCalloc1(2*waterdata->nedge,&edgelist_water);CHKERRQ(ierr); 382 ierr = GetListofEdges_Water(waterdata,edgelist_water);CHKERRQ(ierr); 383 numEdges[1] = waterdata->nedge; 384 numVertices[1] = waterdata->nvertex; 385 } 386 387 /* Get data for the coupling subnetwork */ 388 if (size == 1) { /* TODO: for size > 1, parallel processing coupling is buggy */ 389 nsubnetCouple = 1; 390 numEdgesCouple[0] = 1; 391 392 ierr = PetscMalloc1(4*numEdgesCouple[0],&edgelist_couple);CHKERRQ(ierr); 393 edgelist_couple[0] = 0; edgelist_couple[1] = 4; /* from node: net[0] vertex[4] */ 394 edgelist_couple[2] = 1; edgelist_couple[3] = 0; /* to node: net[1] vertex[0] */ 395 } 396 ierr = PetscLogStagePop();CHKERRQ(ierr); 397 398 /* (2) Create network */ 399 /*--------------------*/ 400 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr); 401 ierr = PetscLogStageRegister("Net Setup",&stage[1]);CHKERRQ(ierr); 402 ierr = PetscLogStagePush(stage[1]);CHKERRQ(ierr); 403 404 ierr = PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);CHKERRQ(ierr); 405 ierr = PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);CHKERRQ(ierr); 406 ierr = PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);CHKERRQ(ierr); 407 408 /* Create an empty network object */ 409 ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr); 410 411 /* Register the components in the network */ 412 ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);CHKERRQ(ierr); 413 ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);CHKERRQ(ierr); 414 ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);CHKERRQ(ierr); 415 ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);CHKERRQ(ierr); 416 417 ierr = DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);CHKERRQ(ierr); 418 ierr = DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);CHKERRQ(ierr); 419 420 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); 421 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 422 423 ierr = DMNetworkSetSizes(networkdm,nsubnet,numVertices,numEdges,nsubnetCouple,numEdgesCouple);CHKERRQ(ierr); 424 425 /* Add edge connectivity */ 426 edgelist[0] = edgelist_power; 427 edgelist[1] = edgelist_water; 428 ierr = DMNetworkSetEdgeList(networkdm,edgelist,&edgelist_couple);CHKERRQ(ierr); 429 430 /* Set up the network layout */ 431 ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr); 432 433 /* Add network components - only process[0] has any data to add */ 434 /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */ 435 genj = 0; loadj = 0; 436 if (rank == 0) { 437 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 438 /* ierr = PetscPrintf(PETSC_COMM_SELF,"Power network: nv %D, ne %D\n",nv,ne);CHKERRQ(ierr); */ 439 for (i = 0; i < ne; i++) { 440 ierr = DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i]);CHKERRQ(ierr); 441 } 442 443 for (i = 0; i < nv; i++) { 444 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i]);CHKERRQ(ierr); 445 if (pfdata->bus[i].ngen) { 446 for (j = 0; j < pfdata->bus[i].ngen; j++) { 447 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++]);CHKERRQ(ierr); 448 } 449 } 450 if (pfdata->bus[i].nload) { 451 for (j=0; j < pfdata->bus[i].nload; j++) { 452 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++]);CHKERRQ(ierr); 453 } 454 } 455 /* Add number of variables */ 456 ierr = DMNetworkAddNumVariables(networkdm,vtx[i],2);CHKERRQ(ierr); 457 } 458 } 459 460 /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */ 461 if (size == 1 || (size > 1 && rank == 1)) { 462 ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 463 /* ierr = PetscPrintf(PETSC_COMM_SELF,"Water network: nv %D, ne %D\n",nv,ne);CHKERRQ(ierr); */ 464 for (i = 0; i < ne; i++) { 465 ierr = DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i]);CHKERRQ(ierr); 466 } 467 468 for (i = 0; i < nv; i++) { 469 ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[i]);CHKERRQ(ierr); 470 /* Add number of variables */ 471 ierr = DMNetworkAddNumVariables(networkdm,vtx[i],1);CHKERRQ(ierr); 472 } 473 } 474 475 /* Set up DM for use */ 476 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 477 if (viewDM) { 478 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");CHKERRQ(ierr); 479 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 480 } 481 482 ierr = DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);CHKERRQ(ierr); 483 if (ne && viewDM) { 484 const PetscInt *cone; 485 PetscInt vid[2]; 486 ierr = DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);CHKERRQ(ierr); 487 488 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);CHKERRQ(ierr); 489 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);CHKERRQ(ierr); 490 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Coupling subnetwork: ne %d; connected vertices %d %d\n",rank,ne,vid[0],vid[1]);CHKERRQ(ierr); 491 } 492 493 /* Free user objects */ 494 if (!rank) { 495 ierr = PetscFree(edgelist_power);CHKERRQ(ierr); 496 ierr = PetscFree(pfdata->bus);CHKERRQ(ierr); 497 ierr = PetscFree(pfdata->gen);CHKERRQ(ierr); 498 ierr = PetscFree(pfdata->branch);CHKERRQ(ierr); 499 ierr = PetscFree(pfdata->load);CHKERRQ(ierr); 500 ierr = PetscFree(pfdata);CHKERRQ(ierr); 501 } 502 if (size == 1 || (size > 1 && rank == 1)) { 503 ierr = PetscFree(edgelist_water);CHKERRQ(ierr); 504 ierr = PetscFree(waterdata->vertex);CHKERRQ(ierr); 505 ierr = PetscFree(waterdata->edge);CHKERRQ(ierr); 506 ierr = PetscFree(waterdata);CHKERRQ(ierr); 507 } 508 if (size == 1) { 509 ierr = PetscFree(edgelist_couple);CHKERRQ(ierr); 510 } 511 512 /* Re-distribute networkdm to multiple processes for better job balance */ 513 if (distribute) { 514 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 515 if (viewDM) { 516 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");CHKERRQ(ierr); 517 ierr = DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 518 } 519 } 520 521 /* Test DMNetworkGetSubnetworkCoupleInfo() */ 522 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr); 523 if (test) { 524 for (i=0; i<nsubnetCouple; i++) { 525 ierr = DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);CHKERRQ(ierr); 526 if (ne && viewDM) { 527 const PetscInt *cone; 528 PetscInt vid[2]; 529 ierr = DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);CHKERRQ(ierr); 530 531 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);CHKERRQ(ierr); 532 ierr = DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);CHKERRQ(ierr); 533 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); 534 } 535 } 536 } 537 538 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 539 ierr = VecDuplicate(X,&F);CHKERRQ(ierr); 540 ierr = DMGetLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 541 542 ierr = PetscLogStagePop();CHKERRQ(ierr); 543 544 /* (3) Setup Solvers */ 545 /*-------------------*/ 546 ierr = PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);CHKERRQ(ierr); 547 ierr = PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);CHKERRQ(ierr); 548 549 ierr = PetscLogStageRegister("SNES Setup",&stage[2]);CHKERRQ(ierr); 550 ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr); 551 552 ierr = SetInitialGuess(networkdm,X,&user);CHKERRQ(ierr); 553 /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 554 555 /* Create coupled snes */ 556 /*-------------------- */ 557 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");CHKERRQ(ierr); 558 user.subsnes_id = nsubnet; 559 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 560 ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr); 561 ierr = SNESSetOptionsPrefix(snes,"coupled_");CHKERRQ(ierr); 562 ierr = SNESSetFunction(snes,F,FormFunction,&user);CHKERRQ(ierr); 563 ierr = SNESMonitorSet(snes,UserMonitor,&user,NULL);CHKERRQ(ierr); 564 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 565 566 if (viewJ) { 567 /* View Jac structure */ 568 ierr = SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);CHKERRQ(ierr); 569 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 570 } 571 572 if (viewX) { 573 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");CHKERRQ(ierr); 574 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 575 } 576 577 if (viewJ) { 578 /* View assembled Jac */ 579 ierr = MatView(Jac,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 580 } 581 582 /* Create snes_power */ 583 /*-------------------*/ 584 ierr = PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");CHKERRQ(ierr); 585 ierr = SetInitialGuess(networkdm,X,&user);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 ierr = SetInitialGuess(networkdm,X,&user);CHKERRQ(ierr); 614 615 user.subsnes_id = 1; 616 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_water);CHKERRQ(ierr); 617 ierr = SNESSetDM(snes_water,networkdm);CHKERRQ(ierr); 618 ierr = SNESSetOptionsPrefix(snes_water,"water_");CHKERRQ(ierr); 619 ierr = SNESSetFunction(snes_water,F,FormFunction,&user);CHKERRQ(ierr); 620 ierr = SNESMonitorSet(snes_water,UserMonitor,&user,NULL);CHKERRQ(ierr); 621 ierr = SNESSetFromOptions(snes_water);CHKERRQ(ierr); 622 623 if (viewX) { 624 ierr = PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");CHKERRQ(ierr); 625 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 626 } 627 ierr = PetscLogStagePop();CHKERRQ(ierr); 628 629 /* (4) Solve */ 630 /*-----------*/ 631 ierr = PetscLogStageRegister("SNES Solve",&stage[3]);CHKERRQ(ierr); 632 ierr = PetscLogStagePush(stage[3]);CHKERRQ(ierr); 633 user.it = 0; 634 reason = SNES_DIVERGED_DTOL; 635 while (user.it < it_max && (PetscInt)reason<0) { 636 #if 0 637 user.subsnes_id = 0; 638 ierr = SNESSolve(snes_power,NULL,X);CHKERRQ(ierr); 639 640 user.subsnes_id = 1; 641 ierr = SNESSolve(snes_water,NULL,X);CHKERRQ(ierr); 642 #endif 643 user.subsnes_id = nsubnet; 644 ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); 645 646 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 647 user.it++; 648 } 649 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);CHKERRQ(ierr); 650 if (viewX) { 651 ierr = PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");CHKERRQ(ierr); 652 ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 653 } 654 ierr = PetscLogStagePop();CHKERRQ(ierr); 655 656 /* Free objects */ 657 /* -------------*/ 658 ierr = VecDestroy(&X);CHKERRQ(ierr); 659 ierr = VecDestroy(&F);CHKERRQ(ierr); 660 ierr = DMRestoreLocalVector(networkdm,&user.localXold);CHKERRQ(ierr); 661 662 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 663 ierr = MatDestroy(&Jac);CHKERRQ(ierr); 664 ierr = SNESDestroy(&snes_power);CHKERRQ(ierr); 665 ierr = SNESDestroy(&snes_water);CHKERRQ(ierr); 666 667 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 668 ierr = PetscFinalize(); 669 return ierr; 670 } 671 672 /*TEST 673 674 build: 675 requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED) 676 depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c 677 678 test: 679 args: -coupled_snes_converged_reason -options_left no 680 localrunfiles: ex1options power/case9.m water/sample1.inp 681 output_file: output/ex1.out 682 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 683 684 test: 685 suffix: 2 686 nsize: 3 687 args: -coupled_snes_converged_reason -options_left no 688 localrunfiles: ex1options power/case9.m water/sample1.inp 689 output_file: output/ex1_2.out 690 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 691 692 test: 693 suffix: 3 694 nsize: 3 695 args: -coupled_snes_converged_reason -options_left no -distribute false 696 localrunfiles: ex1options power/case9.m water/sample1.inp 697 output_file: output/ex1_2.out 698 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 699 700 TEST*/ 701