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