1 2 static char help[] = "This example uses the same problem set up of ex9busdmnetwork.c. \n\ 3 It demonstrates setting and accessing of variables for individual components, instead of \n\ 4 the network vertices (as used in ex9busdmnetwork.c). This is especially useful where vertices \n\ 5 /edges have multiple-components associated with them and one or more components has physics \n\ 6 associated with it. \n\ 7 Input parameters include:\n\ 8 -nc : number of copies of the base case\n\n"; 9 10 /* T 11 Concepts: DMNetwork 12 Concepts: PETSc TS solver 13 14 This example was modified from ex9busdmnetwork.c. 15 */ 16 17 #include <petscts.h> 18 #include <petscdmnetwork.h> 19 20 #define FREQ 60 21 #define W_S (2*PETSC_PI*FREQ) 22 #define NGEN 3 /* No. of generators in the 9 bus system */ 23 #define NLOAD 3 /* No. of loads in the 9 bus system */ 24 #define NBUS 9 /* No. of buses in the 9 bus system */ 25 #define NBRANCH 9 /* No. of branches in the 9 bus system */ 26 27 typedef struct { 28 PetscInt id; /* Bus Number or extended bus name*/ 29 PetscScalar mbase; /* MVA base of the machine */ 30 PetscScalar PG; /* Generator active power output */ 31 PetscScalar QG; /* Generator reactive power output */ 32 33 /* Generator constants */ 34 PetscScalar H; /* Inertia constant */ 35 PetscScalar Rs; /* Stator Resistance */ 36 PetscScalar Xd; /* d-axis reactance */ 37 PetscScalar Xdp; /* d-axis transient reactance */ 38 PetscScalar Xq; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */ 39 PetscScalar Xqp; /* q-axis transient reactance */ 40 PetscScalar Td0p; /* d-axis open circuit time constant */ 41 PetscScalar Tq0p; /* q-axis open circuit time constant */ 42 PetscScalar M; /* M = 2*H/W_S */ 43 PetscScalar D; /* D = 0.1*M */ 44 PetscScalar TM; /* Mechanical Torque */ 45 } Gen; 46 47 typedef struct { 48 /* Exciter system constants */ 49 PetscScalar KA ; /* Voltage regulator gain constant */ 50 PetscScalar TA; /* Voltage regulator time constant */ 51 PetscScalar KE; /* Exciter gain constant */ 52 PetscScalar TE; /* Exciter time constant */ 53 PetscScalar KF; /* Feedback stabilizer gain constant */ 54 PetscScalar TF; /* Feedback stabilizer time constant */ 55 PetscScalar k1,k2; /* calculating the saturation function SE = k1*exp(k2*Efd) */ 56 PetscScalar Vref; /* Voltage regulator voltage setpoint */ 57 } Exc; 58 59 typedef struct { 60 PetscInt id; /* node id */ 61 PetscInt nofgen; /* Number of generators at the bus*/ 62 PetscInt nofload; /* Number of load at the bus*/ 63 PetscScalar yff[2]; /* yff[0]= imaginary part of admittance, yff[1]=real part of admittance*/ 64 PetscScalar vr; /* Real component of bus voltage */ 65 PetscScalar vi; /* Imaginary component of bus voltage */ 66 } Bus; 67 68 /* Load constants 69 We use a composite load model that describes the load and reactive powers at each time instant as follows 70 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i 71 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i 72 where 73 id - index of the load 74 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads 75 ld_alphap,ld_alphap - Percentage contribution (weights) or loads 76 P_D0 - Real power load 77 Q_D0 - Reactive power load 78 Vm(t) - Voltage magnitude at time t 79 Vm0 - Voltage magnitude at t = 0 80 ld_betap, ld_betaq - exponents describing the load model for real and reactive part 81 82 Note: All loads have the same characteristic currently. 83 */ 84 typedef struct { 85 PetscInt id; /* bus id */ 86 PetscInt ld_nsegsp,ld_nsegsq; 87 PetscScalar PD0, QD0; 88 PetscScalar ld_alphap[3]; /* ld_alphap=[1,0,0], an array, not a value, so use *ld_alphap; */ 89 PetscScalar ld_betap[3],ld_alphaq[3],ld_betaq[3]; 90 } Load; 91 92 typedef struct { 93 PetscInt id; /* node id */ 94 PetscScalar yft[2]; /* yft[0]= imaginary part of admittance, yft[1]=real part of admittance*/ 95 } Branch; 96 97 typedef struct { 98 PetscReal tfaulton,tfaultoff; /* Fault on and off times */ 99 PetscReal t; 100 PetscReal t0,tmax; /* initial time and final time */ 101 PetscInt faultbus; /* Fault bus */ 102 PetscScalar Rfault; /* Fault resistance (pu) */ 103 PetscScalar *ybusfault; 104 PetscBool alg_flg; 105 } Userctx; 106 107 /* Used to read data into the DMNetwork components */ 108 PetscErrorCode read_data(PetscInt nc, Gen **pgen,Exc **pexc, Load **pload,Bus **pbus, Branch **pbranch, PetscInt **pedgelist) 109 { 110 PetscErrorCode ierr; 111 PetscInt i,j,row[1],col[2]; 112 PetscInt *edgelist; 113 PetscInt nofgen[9] = {1,1,1,0,0,0,0,0,0}; /* Buses at which generators are incident */ 114 PetscInt nofload[9] = {0,0,0,0,1,1,0,1,0}; /* Buses at which loads are incident */ 115 const PetscScalar *varr; 116 PetscScalar M[3],D[3]; 117 Bus *bus; 118 Branch *branch; 119 Gen *gen; 120 Exc *exc; 121 Load *load; 122 Mat Ybus; 123 Vec V0; 124 125 /*10 parameters*/ 126 /* Generator real and reactive powers (found via loadflow) */ 127 static const PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000}; 128 static const PetscScalar QG[3] = {0.270702180178785,0.066120127797275,-0.108402221791588}; 129 130 /* Generator constants */ 131 static const PetscScalar H[3] = {23.64,6.4,3.01}; /* Inertia constant */ 132 static const PetscScalar Rs[3] = {0.0,0.0,0.0}; /* Stator Resistance */ 133 static const PetscScalar Xd[3] = {0.146,0.8958,1.3125}; /* d-axis reactance */ 134 static const PetscScalar Xdp[3] = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */ 135 static const PetscScalar Xq[3] = {0.4360,0.8645,1.2578}; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */ 136 static const PetscScalar Xqp[3] = {0.0969,0.1969,0.25}; /* q-axis transient reactance */ 137 static const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */ 138 static const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */ 139 140 /* Exciter system constants (8 parameters)*/ 141 static const PetscScalar KA[3] = {20.0,20.0,20.0}; /* Voltage regulartor gain constant */ 142 static const PetscScalar TA[3] = {0.2,0.2,0.2}; /* Voltage regulator time constant */ 143 static const PetscScalar KE[3] = {1.0,1.0,1.0}; /* Exciter gain constant */ 144 static const PetscScalar TE[3] = {0.314,0.314,0.314}; /* Exciter time constant */ 145 static const PetscScalar KF[3] = {0.063,0.063,0.063}; /* Feedback stabilizer gain constant */ 146 static const PetscScalar TF[3] = {0.35,0.35,0.35}; /* Feedback stabilizer time constant */ 147 static const PetscScalar k1[3] = {0.0039,0.0039,0.0039}; 148 static const PetscScalar k2[3] = {1.555,1.555,1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */ 149 150 /* Load constants */ 151 static const PetscScalar PD0[3] = {1.25,0.9,1.0}; 152 static const PetscScalar QD0[3] = {0.5,0.3,0.35}; 153 static const PetscScalar ld_alphaq[3] = {1,0,0}; 154 static const PetscScalar ld_betaq[3] = {2,1,0}; 155 static const PetscScalar ld_betap[3] = {2,1,0}; 156 static const PetscScalar ld_alphap[3] = {1,0,0}; 157 PetscInt ld_nsegsp[3] = {3,3,3}; 158 PetscInt ld_nsegsq[3] = {3,3,3}; 159 PetscViewer Xview,Ybusview; 160 PetscInt neqs_net,m,n; 161 162 PetscFunctionBeginUser; 163 /* Read V0 and Ybus from files */ 164 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"X.bin",FILE_MODE_READ,&Xview);CHKERRQ(ierr); 165 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"Ybus.bin",FILE_MODE_READ,&Ybusview);CHKERRQ(ierr); 166 ierr = VecCreate(PETSC_COMM_SELF,&V0);CHKERRQ(ierr); 167 ierr = VecLoad(V0,Xview);CHKERRQ(ierr); 168 169 ierr = MatCreate(PETSC_COMM_SELF,&Ybus);CHKERRQ(ierr); 170 ierr = MatSetType(Ybus,MATBAIJ);CHKERRQ(ierr); 171 ierr = MatLoad(Ybus,Ybusview);CHKERRQ(ierr); 172 173 /* Destroy unnecessary stuff */ 174 ierr = PetscViewerDestroy(&Xview);CHKERRQ(ierr); 175 ierr = PetscViewerDestroy(&Ybusview);CHKERRQ(ierr); 176 177 ierr = MatGetLocalSize(Ybus,&m,&n);CHKERRQ(ierr); 178 neqs_net = 2*NBUS; /* # eqs. for network subsystem */ 179 if (m != neqs_net || n != neqs_net) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix Ybus is in wrong sizes"); 180 181 M[0] = 2*H[0]/W_S; 182 M[1] = 2*H[1]/W_S; 183 M[2] = 2*H[2]/W_S; 184 D[0] = 0.1*M[0]; 185 D[1] = 0.1*M[1]; 186 D[2] = 0.1*M[2]; 187 188 /* Alocate memory for bus, generators, exciter, loads and branches */ 189 ierr = PetscCalloc5(NBUS*nc,&bus,NGEN*nc,&gen,NLOAD*nc,&load,NBRANCH*nc+(nc-1),&branch,NGEN*nc,&exc);CHKERRQ(ierr); 190 191 ierr = VecGetArrayRead(V0,&varr);CHKERRQ(ierr); 192 193 /* read bus data */ 194 for (i = 0; i < nc; i++) { 195 for (j = 0; j < NBUS; j++) { 196 bus[i*9+j].id = i*9+j; 197 bus[i*9+j].nofgen = nofgen[j]; 198 bus[i*9+j].nofload = nofload[j]; 199 bus[i*9+j].vr = varr[2*j]; 200 bus[i*9+j].vi = varr[2*j+1]; 201 row[0] = 2*j; 202 col[0] = 2*j; 203 col[1] = 2*j+1; 204 /* real and imaginary part of admittance from Ybus into yff */ 205 ierr = MatGetValues(Ybus,1,row,2,col,bus[i*9+j].yff);CHKERRQ(ierr); 206 } 207 } 208 209 /* read generator data */ 210 for (i = 0; i<nc; i++) { 211 for (j = 0; j < NGEN; j++) { 212 gen[i*3+j].id = i*3+j; 213 gen[i*3+j].PG = PG[j]; 214 gen[i*3+j].QG = QG[j]; 215 gen[i*3+j].H = H[j]; 216 gen[i*3+j].Rs = Rs[j]; 217 gen[i*3+j].Xd = Xd[j]; 218 gen[i*3+j].Xdp = Xdp[j]; 219 gen[i*3+j].Xq = Xq[j]; 220 gen[i*3+j].Xqp = Xqp[j]; 221 gen[i*3+j].Td0p = Td0p[j]; 222 gen[i*3+j].Tq0p = Tq0p[j]; 223 gen[i*3+j].M = M[j]; 224 gen[i*3+j].D = D[j]; 225 } 226 } 227 228 for (i = 0; i < nc; i++) { 229 for (j = 0; j < NGEN; j++) { 230 /* exciter system */ 231 exc[i*3+j].KA = KA[j]; 232 exc[i*3+j].TA = TA[j]; 233 exc[i*3+j].KE = KE[j]; 234 exc[i*3+j].TE = TE[j]; 235 exc[i*3+j].KF = KF[j]; 236 exc[i*3+j].TF = TF[j]; 237 exc[i*3+j].k1 = k1[j]; 238 exc[i*3+j].k2 = k2[j]; 239 } 240 } 241 242 /* read load data */ 243 for (i = 0; i<nc; i++) { 244 for (j = 0; j < NLOAD; j++) { 245 load[i*3+j].id = i*3+j; 246 load[i*3+j].PD0 = PD0[j]; 247 load[i*3+j].QD0 = QD0[j]; 248 load[i*3+j].ld_nsegsp = ld_nsegsp[j]; 249 250 load[i*3+j].ld_alphap[0] = ld_alphap[0]; 251 load[i*3+j].ld_alphap[1] = ld_alphap[1]; 252 load[i*3+j].ld_alphap[2] = ld_alphap[2]; 253 254 load[i*3+j].ld_alphaq[0] = ld_alphaq[0]; 255 load[i*3+j].ld_alphaq[1] = ld_alphaq[1]; 256 load[i*3+j].ld_alphaq[2] = ld_alphaq[2]; 257 258 load[i*3+j].ld_betap[0] = ld_betap[0]; 259 load[i*3+j].ld_betap[1] = ld_betap[1]; 260 load[i*3+j].ld_betap[2] = ld_betap[2]; 261 load[i*3+j].ld_nsegsq = ld_nsegsq[j]; 262 263 load[i*3+j].ld_betaq[0] = ld_betaq[0]; 264 load[i*3+j].ld_betaq[1] = ld_betaq[1]; 265 load[i*3+j].ld_betaq[2] = ld_betaq[2]; 266 } 267 } 268 ierr = PetscCalloc1(2*NBRANCH*nc+2*(nc-1),&edgelist);CHKERRQ(ierr); 269 270 /* read edgelist */ 271 for (i = 0; i<nc; i++) { 272 for (j = 0; j < NBRANCH; j++) { 273 switch (j) { 274 case 0: 275 edgelist[i*18+2*j] = 0+9*i; 276 edgelist[i*18+2*j+1] = 3+9*i; 277 break; 278 case 1: 279 edgelist[i*18+2*j] = 1+9*i; 280 edgelist[i*18+2*j+1] = 6+9*i; 281 break; 282 case 2: 283 edgelist[i*18+2*j] = 2+9*i; 284 edgelist[i*18+2*j+1] = 8+9*i; 285 break; 286 case 3: 287 edgelist[i*18+2*j] = 3+9*i; 288 edgelist[i*18+2*j+1] = 4+9*i; 289 break; 290 case 4: 291 edgelist[i*18+2*j] = 3+9*i; 292 edgelist[i*18+2*j+1] = 5+9*i; 293 break; 294 case 5: 295 edgelist[i*18+2*j] = 4+9*i; 296 edgelist[i*18+2*j+1] = 6+9*i; 297 break; 298 case 6: 299 edgelist[i*18+2*j] = 5+9*i; 300 edgelist[i*18+2*j+1] = 8+9*i; 301 break; 302 case 7: 303 edgelist[i*18+2*j] = 6+9*i; 304 edgelist[i*18+2*j+1] = 7+9*i; 305 break; 306 case 8: 307 edgelist[i*18+2*j] = 7+9*i; 308 edgelist[i*18+2*j+1] = 8+9*i; 309 break; 310 default: 311 break; 312 } 313 } 314 } 315 316 /* for connecting last bus of previous network(9*i-1) to first bus of next network(9*i), the branch admittance=-0.0301407+j17.3611 */ 317 for (i = 1; i<nc; i++) { 318 edgelist[18*nc+2*(i-1)] = 8+(i-1)*9; 319 edgelist[18*nc+2*(i-1)+1] = 9*i; 320 321 /* adding admittances to the off-diagonal elements */ 322 branch[9*nc+(i-1)].id = 9*nc+(i-1); 323 branch[9*nc+(i-1)].yft[0] = 17.3611; 324 branch[9*nc+(i-1)].yft[1] = -0.0301407; 325 326 /* subtracting admittances from the diagonal elements */ 327 bus[9*i-1].yff[0] -= 17.3611; 328 bus[9*i-1].yff[1] -= -0.0301407; 329 330 bus[9*i].yff[0] -= 17.3611; 331 bus[9*i].yff[1] -= -0.0301407; 332 } 333 334 /* read branch data */ 335 for (i = 0; i<nc; i++) { 336 for (j = 0; j < NBRANCH; j++) { 337 branch[i*9+j].id = i*9+j; 338 339 row[0] = edgelist[2*j]*2; 340 col[0] = edgelist[2*j+1]*2; 341 col[1] = edgelist[2*j+1]*2+1; 342 ierr = MatGetValues(Ybus,1,row,2,col,branch[i*9+j].yft);CHKERRQ(ierr);/*imaginary part of admittance*/ 343 } 344 } 345 346 *pgen = gen; 347 *pexc = exc; 348 *pload = load; 349 *pbus = bus; 350 *pbranch = branch; 351 *pedgelist = edgelist; 352 353 ierr = VecRestoreArrayRead(V0,&varr);CHKERRQ(ierr); 354 355 /* Destroy unnecessary stuff */ 356 ierr = MatDestroy(&Ybus);CHKERRQ(ierr); 357 ierr = VecDestroy(&V0);CHKERRQ(ierr); 358 PetscFunctionReturn(0); 359 } 360 361 PetscErrorCode SetInitialGuess(DM networkdm, Vec X) 362 { 363 PetscErrorCode ierr; 364 Bus *bus; 365 Gen *gen; 366 Exc *exc; 367 PetscInt v,vStart,vEnd,offset; 368 PetscInt key,numComps,j; 369 PetscBool ghostvtex; 370 Vec localX; 371 PetscScalar *xarr; 372 PetscScalar Vr=0,Vi=0,Vm=0,Vm2; /* Terminal voltage variables */ 373 PetscScalar IGr, IGi; /* Generator real and imaginary current */ 374 PetscScalar Eqp,Edp,delta; /* Generator variables */ 375 PetscScalar Efd=0,RF,VR; /* Exciter variables */ 376 PetscScalar Vd,Vq; /* Generator dq axis voltages */ 377 PetscScalar Id,Iq; /* Generator dq axis currents */ 378 PetscScalar theta; /* Generator phase angle */ 379 PetscScalar SE; 380 void* component; 381 382 PetscFunctionBegin; 383 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 384 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 385 386 ierr = VecSet(X,0.0);CHKERRQ(ierr); 387 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 388 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 389 390 ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr); 391 392 for (v = vStart; v < vEnd; v++) { 393 ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr); 394 if (ghostvtex) continue; 395 396 ierr = DMNetworkGetNumComponents(networkdm,v,&numComps);CHKERRQ(ierr); 397 for (j=0; j < numComps; j++) { 398 ierr = DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL);CHKERRQ(ierr); 399 if (key == 1) { 400 bus = (Bus*)(component); 401 402 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offset);CHKERRQ(ierr); 403 xarr[offset] = bus->vr; 404 xarr[offset+1] = bus->vi; 405 406 Vr = bus->vr; 407 Vi = bus->vi; 408 } else if (key == 2) { 409 gen = (Gen*)(component); 410 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offset);CHKERRQ(ierr); 411 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 412 Vm2 = Vm*Vm; 413 /* Real part of gen current */ 414 IGr = (Vr*gen->PG + Vi*gen->QG)/Vm2; 415 /* Imaginary part of gen current */ 416 IGi = (Vi*gen->PG - Vr*gen->QG)/Vm2; 417 418 /* Machine angle */ 419 delta = atan2(Vi+gen->Xq*IGr,Vr-gen->Xq*IGi); 420 theta = PETSC_PI/2.0 - delta; 421 422 /* d-axis stator current */ 423 Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); 424 425 /* q-axis stator current */ 426 Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); 427 428 Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta); 429 Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta); 430 431 /* d-axis transient EMF */ 432 Edp = Vd + gen->Rs*Id - gen->Xqp*Iq; 433 434 /* q-axis transient EMF */ 435 Eqp = Vq + gen->Rs*Iq + gen->Xdp*Id; 436 437 gen->TM = gen->PG; 438 439 xarr[offset] = Eqp; 440 xarr[offset+1] = Edp; 441 xarr[offset+2] = delta; 442 xarr[offset+3] = W_S; 443 xarr[offset+4] = Id; 444 xarr[offset+5] = Iq; 445 446 Efd = Eqp + (gen->Xd - gen->Xdp)*Id; 447 448 } else if (key == 3) { 449 exc = (Exc*)(component); 450 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offset);CHKERRQ(ierr); 451 452 SE = exc->k1*PetscExpScalar(exc->k2*Efd); 453 VR = exc->KE*Efd + SE; 454 RF = exc->KF*Efd/exc->TF; 455 456 xarr[offset] = Efd; 457 xarr[offset+1] = RF; 458 xarr[offset+2] = VR; 459 460 exc->Vref = Vm + (VR/exc->KA); 461 } 462 } 463 } 464 ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr); 465 ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 466 ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 467 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 468 PetscFunctionReturn(0); 469 } 470 471 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */ 472 PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr,PetscScalar *Fi) 473 { 474 PetscFunctionBegin; 475 *Fr = Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta); 476 *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta); 477 PetscFunctionReturn(0); 478 } 479 480 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */ 481 PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd,PetscScalar *Fq) 482 { 483 PetscFunctionBegin; 484 *Fd = Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta); 485 *Fq = Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta); 486 PetscFunctionReturn(0); 487 } 488 489 /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */ 490 PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx *user) 491 { 492 PetscErrorCode ierr; 493 DM networkdm; 494 Vec localX,localXdot,localF; 495 PetscInt vfrom,vto,offsetfrom,offsetto; 496 PetscInt v,vStart,vEnd,e; 497 PetscScalar *farr; 498 PetscScalar Vd,Vq,SE; 499 const PetscScalar *xarr,*xdotarr; 500 void* component; 501 PetscScalar Vr=0, Vi=0; 502 503 PetscFunctionBegin; 504 ierr = VecSet(F,0.0);CHKERRQ(ierr); 505 506 ierr = TSGetDM(ts,&networkdm);CHKERRQ(ierr); 507 ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr); 508 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 509 ierr = DMGetLocalVector(networkdm,&localXdot);CHKERRQ(ierr); 510 ierr = VecSet(localF,0.0);CHKERRQ(ierr); 511 512 /* update ghost values of localX and localXdot */ 513 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 514 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 515 516 ierr = DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr); 517 ierr = DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr); 518 519 ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); 520 ierr = VecGetArrayRead(localXdot,&xdotarr);CHKERRQ(ierr); 521 ierr = VecGetArray(localF,&farr);CHKERRQ(ierr); 522 523 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 524 525 for (v=vStart; v < vEnd; v++) { 526 PetscInt i,j,offsetbus,offsetgen,offsetexc,key; 527 Bus *bus; 528 Gen *gen; 529 Exc *exc; 530 Load *load; 531 PetscBool ghostvtex; 532 PetscInt numComps; 533 PetscScalar Yffr,Yffi; /* Real and imaginary fault admittances */ 534 PetscScalar Vm,Vm2,Vm0; 535 PetscScalar Vr0=0,Vi0=0; 536 PetscScalar PD,QD; 537 538 ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr); 539 ierr = DMNetworkGetNumComponents(networkdm,v,&numComps);CHKERRQ(ierr); 540 541 for (j = 0; j < numComps; j++) { 542 ierr = DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL);CHKERRQ(ierr); 543 if (key == 1) { 544 PetscInt nconnedges; 545 const PetscInt *connedges; 546 547 bus = (Bus*)(component); 548 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetbus);CHKERRQ(ierr); 549 if (!ghostvtex) { 550 Vr = xarr[offsetbus]; 551 Vi = xarr[offsetbus+1]; 552 553 Yffr = bus->yff[1]; 554 Yffi = bus->yff[0]; 555 556 if (user->alg_flg) { 557 Yffr += user->ybusfault[bus->id*2+1]; 558 Yffi += user->ybusfault[bus->id*2]; 559 } 560 Vr0 = bus->vr; 561 Vi0 = bus->vi; 562 563 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here. 564 The generator current injection, IG, and load current injection, ID are added later 565 */ 566 farr[offsetbus] += Yffi*Vr + Yffr*Vi; /* imaginary current due to diagonal elements */ 567 farr[offsetbus+1] += Yffr*Vr - Yffi*Vi; /* real current due to diagonal elements */ 568 } 569 570 ierr = DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);CHKERRQ(ierr); 571 572 for (i=0; i < nconnedges; i++) { 573 Branch *branch; 574 PetscInt keye; 575 PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti; 576 const PetscInt *cone; 577 578 e = connedges[i]; 579 ierr = DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);CHKERRQ(ierr); 580 581 Yfti = branch->yft[0]; 582 Yftr = branch->yft[1]; 583 584 ierr = DMNetworkGetConnectedVertices(networkdm,e,&cone);CHKERRQ(ierr); 585 586 vfrom = cone[0]; 587 vto = cone[1]; 588 589 ierr = DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom);CHKERRQ(ierr); 590 ierr = DMNetworkGetLocalVecOffset(networkdm,vto,0,&offsetto);CHKERRQ(ierr); 591 592 /* From bus and to bus real and imaginary voltages */ 593 Vfr = xarr[offsetfrom]; 594 Vfi = xarr[offsetfrom+1]; 595 Vtr = xarr[offsetto]; 596 Vti = xarr[offsetto+1]; 597 598 if (vfrom == v) { 599 farr[offsetfrom] += Yftr*Vti + Yfti*Vtr; 600 farr[offsetfrom+1] += Yftr*Vtr - Yfti*Vti; 601 } else { 602 farr[offsetto] += Yftr*Vfi + Yfti*Vfr; 603 farr[offsetto+1] += Yftr*Vfr - Yfti*Vfi; 604 } 605 } 606 } else if (key == 2) { 607 if (!ghostvtex) { 608 PetscScalar Eqp,Edp,delta,w; /* Generator variables */ 609 PetscScalar Efd; /* Exciter field voltage */ 610 PetscScalar Id,Iq; /* Generator dq axis currents */ 611 PetscScalar IGr,IGi,Zdq_inv[4],det; 612 PetscScalar Xd,Xdp,Td0p,Xq,Xqp,Tq0p,TM,D,M,Rs; /* Generator parameters */ 613 614 gen = (Gen*)(component); 615 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetgen);CHKERRQ(ierr); 616 617 /* Generator state variables */ 618 Eqp = xarr[offsetgen]; 619 Edp = xarr[offsetgen+1]; 620 delta = xarr[offsetgen+2]; 621 w = xarr[offsetgen+3]; 622 Id = xarr[offsetgen+4]; 623 Iq = xarr[offsetgen+5]; 624 625 /* Generator parameters */ 626 Xd = gen->Xd; 627 Xdp = gen->Xdp; 628 Td0p = gen->Td0p; 629 Xq = gen->Xq; 630 Xqp = gen->Xqp; 631 Tq0p = gen->Tq0p; 632 TM = gen->TM; 633 D = gen->D; 634 M = gen->M; 635 Rs = gen->Rs; 636 637 ierr = DMNetworkGetLocalVecOffset(networkdm,v,2,&offsetexc);CHKERRQ(ierr); 638 Efd = xarr[offsetexc]; 639 640 /* Generator differential equations */ 641 farr[offsetgen] = (Eqp + (Xd - Xdp)*Id - Efd)/Td0p + xdotarr[offsetgen]; 642 farr[offsetgen+1] = (Edp - (Xq - Xqp)*Iq)/Tq0p + xdotarr[offsetgen+1]; 643 farr[offsetgen+2] = -w + W_S + xdotarr[offsetgen+2]; 644 farr[offsetgen+3] = (-TM + Edp*Id + Eqp*Iq + (Xqp - Xdp)*Id*Iq + D*(w - W_S))/M + xdotarr[offsetgen+3]; 645 646 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 647 648 /* Algebraic equations for stator currents */ 649 det = Rs*Rs + Xdp*Xqp; 650 651 Zdq_inv[0] = Rs/det; 652 Zdq_inv[1] = Xqp/det; 653 Zdq_inv[2] = -Xdp/det; 654 Zdq_inv[3] = Rs/det; 655 656 farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 657 farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 658 659 ierr = dq2ri(Id,Iq,delta,&IGr,&IGi);CHKERRQ(ierr); 660 661 /* Add generator current injection to network */ 662 farr[offsetbus] -= IGi; 663 farr[offsetbus+1] -= IGr; 664 665 } 666 } else if (key == 3) { 667 if (!ghostvtex) { 668 PetscScalar k1,k2,KE,TE,TF,KA,KF,Vref,TA; /* Generator parameters */ 669 PetscScalar Efd,RF,VR; /* Exciter variables */ 670 671 exc = (Exc*)(component); 672 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetexc);CHKERRQ(ierr); 673 674 Efd = xarr[offsetexc]; 675 RF = xarr[offsetexc+1]; 676 VR = xarr[offsetexc+2]; 677 678 k1 = exc->k1; 679 k2 = exc->k2; 680 KE = exc->KE; 681 TE = exc->TE; 682 TF = exc->TF; 683 KA = exc->KA; 684 KF = exc->KF; 685 Vref = exc->Vref; 686 TA = exc->TA; 687 688 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 689 SE = k1*PetscExpScalar(k2*Efd); 690 691 /* Exciter differential equations */ 692 farr[offsetexc] = (KE*Efd + SE - VR)/TE + xdotarr[offsetexc]; 693 farr[offsetexc+1] = (RF - KF*Efd/TF)/TF + xdotarr[offsetexc+1]; 694 farr[offsetexc+2] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[offsetexc+2]; 695 696 } 697 } else if (key ==4) { 698 if (!ghostvtex) { 699 PetscInt k; 700 PetscInt ld_nsegsp; 701 PetscInt ld_nsegsq; 702 PetscScalar *ld_alphap; 703 PetscScalar *ld_betap,*ld_alphaq,*ld_betaq,PD0, QD0, IDr,IDi; 704 705 load = (Load*)(component); 706 707 /* Load Parameters */ 708 ld_nsegsp = load->ld_nsegsp; 709 ld_alphap = load->ld_alphap; 710 ld_betap = load->ld_betap; 711 ld_nsegsq = load->ld_nsegsq; 712 ld_alphaq = load->ld_alphaq; 713 ld_betaq = load->ld_betaq; 714 PD0 = load->PD0; 715 QD0 = load->QD0; 716 717 Vr = xarr[offsetbus]; /* Real part of generator terminal voltage */ 718 Vi = xarr[offsetbus+1]; /* Imaginary part of the generator terminal voltage */ 719 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 720 Vm2 = Vm*Vm; 721 Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0); 722 PD = QD = 0.0; 723 for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]); 724 for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 725 726 /* Load currents */ 727 IDr = (PD*Vr + QD*Vi)/Vm2; 728 IDi = (-QD*Vr + PD*Vi)/Vm2; 729 730 /* Load current contribution to the network */ 731 farr[offsetbus] += IDi; 732 farr[offsetbus+1] += IDr; 733 } 734 } 735 } 736 } 737 738 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 739 ierr = VecRestoreArrayRead(localXdot,&xdotarr);CHKERRQ(ierr); 740 ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); 741 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 742 ierr = DMRestoreLocalVector(networkdm,&localXdot);CHKERRQ(ierr); 743 744 ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 745 ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 746 ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr); 747 PetscFunctionReturn(0); 748 } 749 750 /* This function is used for solving the algebraic system only during fault on and 751 off times. It computes the entire F and then zeros out the part corresponding to 752 differential equations 753 F = [0;g(y)]; 754 */ 755 PetscErrorCode AlgFunction (SNES snes, Vec X, Vec F, void *ctx) 756 { 757 PetscErrorCode ierr; 758 DM networkdm; 759 Vec localX,localF; 760 PetscInt vfrom,vto,offsetfrom,offsetto; 761 PetscInt v,vStart,vEnd,e; 762 PetscScalar *farr; 763 Userctx *user=(Userctx*)ctx; 764 const PetscScalar *xarr; 765 void* component; 766 PetscScalar Vr=0,Vi=0; 767 768 PetscFunctionBegin; 769 ierr = VecSet(F,0.0);CHKERRQ(ierr); 770 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 771 ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr); 772 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 773 ierr = VecSet(localF,0.0);CHKERRQ(ierr); 774 775 /* update ghost values of locaX and locaXdot */ 776 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 777 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 778 779 ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); 780 ierr = VecGetArray(localF,&farr);CHKERRQ(ierr); 781 782 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 783 784 for (v=vStart; v < vEnd; v++) { 785 PetscInt i,j,offsetbus,offsetgen,key,numComps; 786 PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0=0, Vi0=0, PD, QD; 787 Bus *bus; 788 Gen *gen; 789 Load *load; 790 PetscBool ghostvtex; 791 792 ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr); 793 ierr = DMNetworkGetNumComponents(networkdm,v,&numComps);CHKERRQ(ierr); 794 795 for (j = 0; j < numComps; j++) { 796 ierr = DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL);CHKERRQ(ierr); 797 if (key == 1) { 798 PetscInt nconnedges; 799 const PetscInt *connedges; 800 801 bus = (Bus*)(component); 802 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetbus);CHKERRQ(ierr); 803 if (!ghostvtex) { 804 Vr = xarr[offsetbus]; 805 Vi = xarr[offsetbus+1]; 806 807 Yffr = bus->yff[1]; 808 Yffi = bus->yff[0]; 809 if (user->alg_flg) { 810 Yffr += user->ybusfault[bus->id*2+1]; 811 Yffi += user->ybusfault[bus->id*2]; 812 } 813 Vr0 = bus->vr; 814 Vi0 = bus->vi; 815 816 farr[offsetbus] += Yffi*Vr + Yffr*Vi; 817 farr[offsetbus+1] += Yffr*Vr - Yffi*Vi; 818 } 819 ierr = DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);CHKERRQ(ierr); 820 821 for (i=0; i < nconnedges; i++) { 822 Branch *branch; 823 PetscInt keye; 824 PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti; 825 const PetscInt *cone; 826 827 e = connedges[i]; 828 ierr = DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);CHKERRQ(ierr); 829 830 Yfti = branch->yft[0]; 831 Yftr = branch->yft[1]; 832 833 ierr = DMNetworkGetConnectedVertices(networkdm,e,&cone);CHKERRQ(ierr); 834 vfrom = cone[0]; 835 vto = cone[1]; 836 837 ierr = DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom);CHKERRQ(ierr); 838 ierr = DMNetworkGetLocalVecOffset(networkdm,vto,0,&offsetto);CHKERRQ(ierr); 839 840 /*From bus and to bus real and imaginary voltages */ 841 Vfr = xarr[offsetfrom]; 842 Vfi = xarr[offsetfrom+1]; 843 Vtr = xarr[offsetto]; 844 Vti = xarr[offsetto+1]; 845 846 if (vfrom == v) { 847 farr[offsetfrom] += Yftr*Vti + Yfti*Vtr; 848 farr[offsetfrom+1] += Yftr*Vtr - Yfti*Vti; 849 } else { 850 farr[offsetto] += Yftr*Vfi + Yfti*Vfr; 851 farr[offsetto+1] += Yftr*Vfr - Yfti*Vfi; 852 } 853 } 854 } else if (key == 2) { 855 if (!ghostvtex) { 856 PetscScalar Eqp,Edp,delta; /* Generator variables */ 857 PetscScalar Id,Iq; /* Generator dq axis currents */ 858 PetscScalar Vd,Vq,IGr,IGi,Zdq_inv[4],det; 859 PetscScalar Xdp,Xqp,Rs; /* Generator parameters */ 860 861 gen = (Gen*)(component); 862 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetgen);CHKERRQ(ierr); 863 864 /* Generator state variables */ 865 Eqp = xarr[offsetgen]; 866 Edp = xarr[offsetgen+1]; 867 delta = xarr[offsetgen+2]; 868 /* w = xarr[idx+3]; not being used */ 869 Id = xarr[offsetgen+4]; 870 Iq = xarr[offsetgen+5]; 871 872 /* Generator parameters */ 873 Xdp = gen->Xdp; 874 Xqp = gen->Xqp; 875 Rs = gen->Rs; 876 877 /* Set generator differential equation residual functions to zero */ 878 farr[offsetgen] = 0; 879 farr[offsetgen+1] = 0; 880 farr[offsetgen+2] = 0; 881 farr[offsetgen+3] = 0; 882 883 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 884 885 /* Algebraic equations for stator currents */ 886 det = Rs*Rs + Xdp*Xqp; 887 888 Zdq_inv[0] = Rs/det; 889 Zdq_inv[1] = Xqp/det; 890 Zdq_inv[2] = -Xdp/det; 891 Zdq_inv[3] = Rs/det; 892 893 farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 894 farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 895 896 /* Add generator current injection to network */ 897 ierr = dq2ri(Id,Iq,delta,&IGr,&IGi);CHKERRQ(ierr); 898 899 farr[offsetbus] -= IGi; 900 farr[offsetbus+1] -= IGr; 901 902 /* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */ 903 904 } 905 } else if (key == 3) { 906 if (!ghostvtex) { 907 PetscInt offsetexc; 908 ierr = DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetexc);CHKERRQ(ierr); 909 /* Set exciter differential equation residual functions equal to zero*/ 910 farr[offsetexc] = 0; 911 farr[offsetexc+1] = 0; 912 farr[offsetexc+2] = 0; 913 } 914 } else if (key == 4) { 915 if (!ghostvtex) { 916 PetscInt k,ld_nsegsp,ld_nsegsq; 917 PetscScalar *ld_alphap,*ld_betap,*ld_alphaq,*ld_betaq,PD0,QD0,IDr,IDi; 918 919 load = (Load*)(component); 920 921 /* Load Parameters */ 922 ld_nsegsp = load->ld_nsegsp; 923 ld_alphap = load->ld_alphap; 924 ld_betap = load->ld_betap; 925 ld_nsegsq = load->ld_nsegsq; 926 ld_alphaq = load->ld_alphaq; 927 ld_betaq = load->ld_betaq; 928 929 PD0 = load->PD0; 930 QD0 = load->QD0; 931 932 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 933 Vm2 = Vm*Vm; 934 Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0); 935 PD = QD = 0.0; 936 for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]); 937 for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 938 939 /* Load currents */ 940 IDr = (PD*Vr + QD*Vi)/Vm2; 941 IDi = (-QD*Vr + PD*Vi)/Vm2; 942 943 farr[offsetbus] += IDi; 944 farr[offsetbus+1] += IDr; 945 } 946 } 947 } 948 } 949 950 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 951 ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); 952 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 953 954 ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 955 ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 956 ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr); 957 PetscFunctionReturn(0); 958 } 959 960 int main(int argc,char ** argv) 961 { 962 PetscErrorCode ierr; 963 PetscInt i,j,*edgelist= NULL,eStart,eEnd,vStart,vEnd; 964 PetscInt genj,excj,loadj,componentkey[5]; 965 PetscInt nc = 1; /* No. of copies (default = 1) */ 966 PetscMPIInt size,rank; 967 Vec X,F_alg; 968 TS ts; 969 SNES snes_alg,snes; 970 Bus *bus; 971 Branch *branch; 972 Gen *gen; 973 Exc *exc; 974 Load *load; 975 DM networkdm; 976 #if defined(PETSC_USE_LOG) 977 PetscLogStage stage1; 978 #endif 979 Userctx user; 980 KSP ksp; 981 PC pc; 982 PetscInt numEdges = 0,numVertices = 0; 983 984 ierr = PetscInitialize(&argc,&argv,"ex9busnetworkops",help);if (ierr) return ierr; 985 ierr = PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);CHKERRQ(ierr); 986 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 987 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 988 989 /* Read initial voltage vector and Ybus */ 990 if (!rank) { 991 ierr = read_data(nc,&gen,&exc,&load,&bus,&branch,&edgelist);CHKERRQ(ierr); 992 } 993 994 ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr); 995 ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(Branch),&componentkey[0]);CHKERRQ(ierr); 996 ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(Bus),&componentkey[1]);CHKERRQ(ierr); 997 ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(Gen),&componentkey[2]);CHKERRQ(ierr); 998 ierr = DMNetworkRegisterComponent(networkdm,"excstruct",sizeof(Exc),&componentkey[3]);CHKERRQ(ierr); 999 ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(Load),&componentkey[4]);CHKERRQ(ierr); 1000 1001 ierr = PetscLogStageRegister("Create network",&stage1);CHKERRQ(ierr); 1002 ierr = PetscLogStagePush(stage1);CHKERRQ(ierr); 1003 1004 /* Set local number of nodes and edges and edge connectivity */ 1005 if (!rank) { 1006 numVertices = NBUS*nc; numEdges = NBRANCH*nc+(nc-1); 1007 } 1008 ierr = DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);CHKERRQ(ierr); 1009 ierr = DMNetworkAddSubnetwork(networkdm,NULL,numVertices,numEdges,edgelist,NULL);CHKERRQ(ierr); 1010 1011 /* Set up the network layout */ 1012 ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr); 1013 1014 if (!rank) { 1015 ierr = PetscFree(edgelist);CHKERRQ(ierr); 1016 } 1017 1018 /* Add network components (physical parameters of nodes and branches) and number of variables */ 1019 if (!rank) { 1020 ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr); 1021 genj=0; loadj=0; excj=0; 1022 for (i = eStart; i < eEnd; i++) { 1023 ierr = DMNetworkAddComponent(networkdm,i,componentkey[0],&branch[i-eStart],0);CHKERRQ(ierr); 1024 } 1025 1026 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 1027 1028 for (i = vStart; i < vEnd; i++) { 1029 ierr = DMNetworkAddComponent(networkdm,i,componentkey[1],&bus[i-vStart],2);CHKERRQ(ierr); 1030 if (bus[i-vStart].nofgen) { 1031 for (j = 0; j < bus[i-vStart].nofgen; j++) { 1032 /* Add generator */ 1033 ierr = DMNetworkAddComponent(networkdm,i,componentkey[2],&gen[genj++],6);CHKERRQ(ierr); 1034 /* Add exciter */ 1035 ierr = DMNetworkAddComponent(networkdm,i,componentkey[3],&exc[excj++],3);CHKERRQ(ierr); 1036 } 1037 } 1038 if (bus[i-vStart].nofload) { 1039 for (j=0; j < bus[i-vStart].nofload; j++) { 1040 ierr = DMNetworkAddComponent(networkdm,i,componentkey[4],&load[loadj++],0);CHKERRQ(ierr); 1041 } 1042 } 1043 } 1044 } 1045 1046 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 1047 1048 if (!rank) { 1049 ierr = PetscFree5(bus,gen,load,branch,exc);CHKERRQ(ierr); 1050 } 1051 1052 /* for parallel options: Network partitioning and distribution of data */ 1053 if (size > 1) { 1054 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 1055 } 1056 ierr = PetscLogStagePop();CHKERRQ(ierr); 1057 1058 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 1059 1060 ierr = SetInitialGuess(networkdm,X);CHKERRQ(ierr); 1061 1062 /* Options for fault simulation */ 1063 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");CHKERRQ(ierr); 1064 user.tfaulton = 0.02; 1065 user.tfaultoff = 0.05; 1066 user.Rfault = 0.0001; 1067 user.faultbus = 8; 1068 ierr = PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL);CHKERRQ(ierr); 1069 ierr = PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL);CHKERRQ(ierr); 1070 ierr = PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL);CHKERRQ(ierr); 1071 user.t0 = 0.0; 1072 user.tmax = 0.1; 1073 ierr = PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL);CHKERRQ(ierr); 1074 ierr = PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL);CHKERRQ(ierr); 1075 1076 ierr = PetscMalloc1(18*nc,&user.ybusfault);CHKERRQ(ierr); 1077 for (i = 0; i < 18*nc; i++) { 1078 user.ybusfault[i] = 0; 1079 } 1080 user.ybusfault[user.faultbus*2+1] = 1/user.Rfault; 1081 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1082 1083 /* Setup TS solver */ 1084 /*--------------------------------------------------------*/ 1085 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 1086 ierr = TSSetDM(ts,(DM)networkdm);CHKERRQ(ierr); 1087 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 1088 1089 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1090 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 1091 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 1092 ierr = PCSetType(pc,PCBJACOBI);CHKERRQ(ierr); 1093 1094 ierr = TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&user);CHKERRQ(ierr); 1095 ierr = TSSetMaxTime(ts,user.tfaulton);CHKERRQ(ierr); 1096 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 1097 ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr); 1098 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 1099 1100 /*user.alg_flg = PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix. 1101 eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */ 1102 user.alg_flg = PETSC_FALSE; 1103 1104 /* Prefault period */ 1105 if (!rank) { 1106 ierr = PetscPrintf(PETSC_COMM_SELF,"... (1) Prefault period ... \n");CHKERRQ(ierr); 1107 } 1108 1109 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 1110 ierr = TSSetUp(ts);CHKERRQ(ierr); 1111 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1112 1113 /* Create the nonlinear solver for solving the algebraic system */ 1114 ierr = VecDuplicate(X,&F_alg);CHKERRQ(ierr); 1115 1116 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_alg);CHKERRQ(ierr); 1117 ierr = SNESSetDM(snes_alg,(DM)networkdm);CHKERRQ(ierr); 1118 ierr = SNESSetFunction(snes_alg,F_alg,AlgFunction,&user);CHKERRQ(ierr); 1119 ierr = SNESSetOptionsPrefix(snes_alg,"alg_");CHKERRQ(ierr); 1120 ierr = SNESSetFromOptions(snes_alg);CHKERRQ(ierr); 1121 1122 /* Apply disturbance - resistive fault at user.faultbus */ 1123 /* This is done by adding shunt conductance to the diagonal location 1124 in the Ybus matrix */ 1125 user.alg_flg = PETSC_TRUE; 1126 1127 /* Solve the algebraic equations */ 1128 if (!rank) { 1129 ierr = PetscPrintf(PETSC_COMM_SELF,"\n... (2) Apply disturbance, solve algebraic equations ... \n");CHKERRQ(ierr); 1130 } 1131 ierr = SNESSolve(snes_alg,NULL,X);CHKERRQ(ierr); 1132 1133 /* Disturbance period */ 1134 ierr = TSSetTime(ts,user.tfaulton);CHKERRQ(ierr); 1135 ierr = TSSetMaxTime(ts,user.tfaultoff);CHKERRQ(ierr); 1136 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 1137 ierr = TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&user);CHKERRQ(ierr); 1138 1139 user.alg_flg = PETSC_TRUE; 1140 if (!rank) { 1141 ierr = PetscPrintf(PETSC_COMM_SELF,"\n... (3) Disturbance period ... \n");CHKERRQ(ierr); 1142 } 1143 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1144 1145 /* Remove the fault */ 1146 ierr = SNESSetFunction(snes_alg,F_alg,AlgFunction,&user);CHKERRQ(ierr); 1147 1148 user.alg_flg = PETSC_FALSE; 1149 /* Solve the algebraic equations */ 1150 if (!rank) { 1151 ierr = PetscPrintf(PETSC_COMM_SELF,"\n... (4) Remove fault, solve algebraic equations ... \n");CHKERRQ(ierr); 1152 } 1153 ierr = SNESSolve(snes_alg,NULL,X);CHKERRQ(ierr); 1154 ierr = SNESDestroy(&snes_alg);CHKERRQ(ierr); 1155 1156 /* Post-disturbance period */ 1157 ierr = TSSetTime(ts,user.tfaultoff);CHKERRQ(ierr); 1158 ierr = TSSetMaxTime(ts,user.tmax);CHKERRQ(ierr); 1159 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 1160 ierr = TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&user);CHKERRQ(ierr); 1161 1162 user.alg_flg = PETSC_FALSE; 1163 if (!rank) { 1164 ierr = PetscPrintf(PETSC_COMM_SELF,"\n... (5) Post-disturbance period ... \n");CHKERRQ(ierr); 1165 } 1166 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1167 1168 ierr = PetscFree(user.ybusfault);CHKERRQ(ierr); 1169 ierr = VecDestroy(&F_alg);CHKERRQ(ierr); 1170 ierr = VecDestroy(&X);CHKERRQ(ierr); 1171 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 1172 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1173 ierr = PetscFinalize(); 1174 return ierr; 1175 } 1176 1177 /*TEST 1178 1179 build: 1180 requires: double !complex !define(PETSC_USE_64BIT_INDICES) 1181 1182 test: 1183 args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason 1184 localrunfiles: X.bin Ybus.bin ex9busnetworkops 1185 1186 test: 1187 suffix: 2 1188 nsize: 2 1189 args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason 1190 localrunfiles: X.bin Ybus.bin ex9busnetworkops 1191 1192 TEST*/ 1193