1 2 static char help[] = "Sensitivity analysis applied in power grid stability analysis of WECC 9 bus system.\n\ 3 This example is based on the 9-bus (node) example given in the book Power\n\ 4 Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\ 5 The power grid in this example consists of 9 buses (nodes), 3 generators,\n\ 6 3 loads, and 9 transmission lines. The network equations are written\n\ 7 in current balance form using rectangular coordiantes.\n\n"; 8 9 /* 10 The equations for the stability analysis are described by the DAE. See ex9bus.c for details. 11 The system has 'jumps' due to faults, thus the time interval is split into multiple sections, and TSSolve is called for each of them. But TSAdjointSolve only needs to be called once since the whole trajectory has been saved in the forward run. 12 The code computes the sensitivity of a final state w.r.t. initial conditions. 13 */ 14 15 #include <petscts.h> 16 #include <petscdm.h> 17 #include <petscdmda.h> 18 #include <petscdmcomposite.h> 19 20 #define freq 60 21 #define w_s (2*PETSC_PI*freq) 22 23 /* Sizes and indices */ 24 const PetscInt nbus = 9; /* Number of network buses */ 25 const PetscInt ngen = 3; /* Number of generators */ 26 const PetscInt nload = 3; /* Number of loads */ 27 const PetscInt gbus[3] = {0,1,2}; /* Buses at which generators are incident */ 28 const PetscInt lbus[3] = {4,5,7}; /* Buses at which loads are incident */ 29 30 /* Generator real and reactive powers (found via loadflow) */ 31 const PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000}; 32 const PetscScalar QG[3] = {0.270702180178785,0.066120127797275,-0.108402221791588}; 33 /* Generator constants */ 34 const PetscScalar H[3] = {23.64,6.4,3.01}; /* Inertia constant */ 35 const PetscScalar Rs[3] = {0.0,0.0,0.0}; /* Stator Resistance */ 36 const PetscScalar Xd[3] = {0.146,0.8958,1.3125}; /* d-axis reactance */ 37 const PetscScalar Xdp[3] = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */ 38 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 */ 39 const PetscScalar Xqp[3] = {0.0969,0.1969,0.25}; /* q-axis transient reactance */ 40 const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */ 41 const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */ 42 PetscScalar M[3]; /* M = 2*H/w_s */ 43 PetscScalar D[3]; /* D = 0.1*M */ 44 45 PetscScalar TM[3]; /* Mechanical Torque */ 46 /* Exciter system constants */ 47 const PetscScalar KA[3] = {20.0,20.0,20.0}; /* Voltage regulartor gain constant */ 48 const PetscScalar TA[3] = {0.2,0.2,0.2}; /* Voltage regulator time constant */ 49 const PetscScalar KE[3] = {1.0,1.0,1.0}; /* Exciter gain constant */ 50 const PetscScalar TE[3] = {0.314,0.314,0.314}; /* Exciter time constant */ 51 const PetscScalar KF[3] = {0.063,0.063,0.063}; /* Feedback stabilizer gain constant */ 52 const PetscScalar TF[3] = {0.35,0.35,0.35}; /* Feedback stabilizer time constant */ 53 const PetscScalar k1[3] = {0.0039,0.0039,0.0039}; 54 const PetscScalar k2[3] = {1.555,1.555,1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */ 55 56 PetscScalar Vref[3]; 57 /* Load constants 58 We use a composite load model that describes the load and reactive powers at each time instant as follows 59 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i 60 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i 61 where 62 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads 63 ld_alphap,ld_alphap - Percentage contribution (weights) or loads 64 P_D0 - Real power load 65 Q_D0 - Reactive power load 66 V_m(t) - Voltage magnitude at time t 67 V_m0 - Voltage magnitude at t = 0 68 ld_betap, ld_betaq - exponents describing the load model for real and reactive part 69 70 Note: All loads have the same characteristic currently. 71 */ 72 const PetscScalar PD0[3] = {1.25,0.9,1.0}; 73 const PetscScalar QD0[3] = {0.5,0.3,0.35}; 74 const PetscInt ld_nsegsp[3] = {3,3,3}; 75 const PetscScalar ld_alphap[3] = {1.0,0.0,0.0}; 76 const PetscScalar ld_betap[3] = {2.0,1.0,0.0}; 77 const PetscInt ld_nsegsq[3] = {3,3,3}; 78 const PetscScalar ld_alphaq[3] = {1.0,0.0,0.0}; 79 const PetscScalar ld_betaq[3] = {2.0,1.0,0.0}; 80 81 typedef struct { 82 DM dmgen, dmnet; /* DMs to manage generator and network subsystem */ 83 DM dmpgrid; /* Composite DM to manage the entire power grid */ 84 Mat Ybus; /* Network admittance matrix */ 85 Vec V0; /* Initial voltage vector (Power flow solution) */ 86 PetscReal tfaulton,tfaultoff; /* Fault on and off times */ 87 PetscInt faultbus; /* Fault bus */ 88 PetscScalar Rfault; 89 PetscReal t0,tmax; 90 PetscInt neqs_gen,neqs_net,neqs_pgrid; 91 PetscBool alg_flg; 92 PetscReal t; 93 IS is_diff; /* indices for differential equations */ 94 IS is_alg; /* indices for algebraic equations */ 95 } Userctx; 96 97 98 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */ 99 PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr, PetscScalar *Fi) 100 { 101 PetscFunctionBegin; 102 *Fr = Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta); 103 *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta); 104 PetscFunctionReturn(0); 105 } 106 107 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */ 108 PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd, PetscScalar *Fq) 109 { 110 PetscFunctionBegin; 111 *Fd = Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta); 112 *Fq = Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta); 113 PetscFunctionReturn(0); 114 } 115 116 PetscErrorCode SetInitialGuess(Vec X,Userctx *user) 117 { 118 PetscErrorCode ierr; 119 Vec Xgen,Xnet; 120 PetscScalar *xgen,*xnet; 121 PetscInt i,idx=0; 122 PetscScalar Vr,Vi,IGr,IGi,Vm,Vm2; 123 PetscScalar Eqp,Edp,delta; 124 PetscScalar Efd,RF,VR; /* Exciter variables */ 125 PetscScalar Id,Iq; /* Generator dq axis currents */ 126 PetscScalar theta,Vd,Vq,SE; 127 128 PetscFunctionBegin; 129 M[0] = 2*H[0]/w_s; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s; 130 D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2]; 131 132 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 133 134 /* Network subsystem initialization */ 135 ierr = VecCopy(user->V0,Xnet);CHKERRQ(ierr); 136 137 /* Generator subsystem initialization */ 138 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 139 ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr); 140 141 for (i=0; i < ngen; i++) { 142 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 143 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 144 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 145 IGr = (Vr*PG[i] + Vi*QG[i])/Vm2; 146 IGi = (Vi*PG[i] - Vr*QG[i])/Vm2; 147 148 delta = PetscAtan2Real(Vi+Xq[i]*IGr,Vr-Xq[i]*IGi); /* Machine angle */ 149 150 theta = PETSC_PI/2.0 - delta; 151 152 Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); /* d-axis stator current */ 153 Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); /* q-axis stator current */ 154 155 Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta); 156 Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta); 157 158 Edp = Vd + Rs[i]*Id - Xqp[i]*Iq; /* d-axis transient EMF */ 159 Eqp = Vq + Rs[i]*Iq + Xdp[i]*Id; /* q-axis transient EMF */ 160 161 TM[i] = PG[i]; 162 163 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */ 164 xgen[idx] = Eqp; 165 xgen[idx+1] = Edp; 166 xgen[idx+2] = delta; 167 xgen[idx+3] = w_s; 168 169 idx = idx + 4; 170 171 xgen[idx] = Id; 172 xgen[idx+1] = Iq; 173 174 idx = idx + 2; 175 176 /* Exciter */ 177 Efd = Eqp + (Xd[i] - Xdp[i])*Id; 178 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 179 VR = KE[i]*Efd + SE; 180 RF = KF[i]*Efd/TF[i]; 181 182 xgen[idx] = Efd; 183 xgen[idx+1] = RF; 184 xgen[idx+2] = VR; 185 186 Vref[i] = Vm + (VR/KA[i]); 187 188 idx = idx + 3; 189 } 190 191 ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr); 192 ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr); 193 194 /* ierr = VecView(Xgen,0);CHKERRQ(ierr); */ 195 ierr = DMCompositeGather(user->dmpgrid,INSERT_VALUES,X,Xgen,Xnet);CHKERRQ(ierr); 196 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 /* Computes F = [f(x,y);g(x,y)] */ 201 PetscErrorCode ResidualFunction(SNES snes,Vec X, Vec F, Userctx *user) 202 { 203 PetscErrorCode ierr; 204 Vec Xgen,Xnet,Fgen,Fnet; 205 PetscScalar *xgen,*xnet,*fgen,*fnet; 206 PetscInt i,idx=0; 207 PetscScalar Vr,Vi,Vm,Vm2; 208 PetscScalar Eqp,Edp,delta,w; /* Generator variables */ 209 PetscScalar Efd,RF,VR; /* Exciter variables */ 210 PetscScalar Id,Iq; /* Generator dq axis currents */ 211 PetscScalar Vd,Vq,SE; 212 PetscScalar IGr,IGi,IDr,IDi; 213 PetscScalar Zdq_inv[4],det; 214 PetscScalar PD,QD,Vm0,*v0; 215 PetscInt k; 216 217 PetscFunctionBegin; 218 ierr = VecZeroEntries(F);CHKERRQ(ierr); 219 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 220 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Fgen,&Fnet);CHKERRQ(ierr); 221 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 222 ierr = DMCompositeScatter(user->dmpgrid,F,Fgen,Fnet);CHKERRQ(ierr); 223 224 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here. 225 The generator current injection, IG, and load current injection, ID are added later 226 */ 227 /* Note that the values in Ybus are stored assuming the imaginary current balance 228 equation is ordered first followed by real current balance equation for each bus. 229 Thus imaginary current contribution goes in location 2*i, and 230 real current contribution in 2*i+1 231 */ 232 ierr = MatMult(user->Ybus,Xnet,Fnet);CHKERRQ(ierr); 233 234 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 235 ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr); 236 ierr = VecGetArray(Fgen,&fgen);CHKERRQ(ierr); 237 ierr = VecGetArray(Fnet,&fnet);CHKERRQ(ierr); 238 239 /* Generator subsystem */ 240 for (i=0; i < ngen; i++) { 241 Eqp = xgen[idx]; 242 Edp = xgen[idx+1]; 243 delta = xgen[idx+2]; 244 w = xgen[idx+3]; 245 Id = xgen[idx+4]; 246 Iq = xgen[idx+5]; 247 Efd = xgen[idx+6]; 248 RF = xgen[idx+7]; 249 VR = xgen[idx+8]; 250 251 /* Generator differential equations */ 252 fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; 253 fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; 254 fgen[idx+2] = -w + w_s; 255 fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; 256 257 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 258 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 259 260 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 261 /* Algebraic equations for stator currents */ 262 263 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 264 265 Zdq_inv[0] = Rs[i]/det; 266 Zdq_inv[1] = Xqp[i]/det; 267 Zdq_inv[2] = -Xdp[i]/det; 268 Zdq_inv[3] = Rs[i]/det; 269 270 fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 271 fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 272 273 /* Add generator current injection to network */ 274 ierr = dq2ri(Id,Iq,delta,&IGr,&IGi);CHKERRQ(ierr); 275 276 fnet[2*gbus[i]] -= IGi; 277 fnet[2*gbus[i]+1] -= IGr; 278 279 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 280 281 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 282 283 /* Exciter differential equations */ 284 fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; 285 fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; 286 fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; 287 288 idx = idx + 9; 289 } 290 291 ierr = VecGetArray(user->V0,&v0);CHKERRQ(ierr); 292 for (i=0; i < nload; i++) { 293 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 294 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 295 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 296 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 297 PD = QD = 0.0; 298 for (k=0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 299 for (k=0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 300 301 /* Load currents */ 302 IDr = (PD*Vr + QD*Vi)/Vm2; 303 IDi = (-QD*Vr + PD*Vi)/Vm2; 304 305 fnet[2*lbus[i]] += IDi; 306 fnet[2*lbus[i]+1] += IDr; 307 } 308 ierr = VecRestoreArray(user->V0,&v0);CHKERRQ(ierr); 309 310 ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr); 311 ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr); 312 ierr = VecRestoreArray(Fgen,&fgen);CHKERRQ(ierr); 313 ierr = VecRestoreArray(Fnet,&fnet);CHKERRQ(ierr); 314 315 ierr = DMCompositeGather(user->dmpgrid,INSERT_VALUES,F,Fgen,Fnet);CHKERRQ(ierr); 316 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 317 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Fgen,&Fnet);CHKERRQ(ierr); 318 PetscFunctionReturn(0); 319 } 320 321 /* \dot{x} - f(x,y) 322 g(x,y) = 0 323 */ 324 PetscErrorCode IFunction(TS ts,PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user) 325 { 326 PetscErrorCode ierr; 327 SNES snes; 328 PetscScalar *f; 329 const PetscScalar *xdot; 330 PetscInt i; 331 332 PetscFunctionBegin; 333 user->t = t; 334 335 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 336 ierr = ResidualFunction(snes,X,F,user);CHKERRQ(ierr); 337 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 338 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 339 for (i=0;i < ngen;i++) { 340 f[9*i] += xdot[9*i]; 341 f[9*i+1] += xdot[9*i+1]; 342 f[9*i+2] += xdot[9*i+2]; 343 f[9*i+3] += xdot[9*i+3]; 344 f[9*i+6] += xdot[9*i+6]; 345 f[9*i+7] += xdot[9*i+7]; 346 f[9*i+8] += xdot[9*i+8]; 347 } 348 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 349 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 350 PetscFunctionReturn(0); 351 } 352 353 /* This function is used for solving the algebraic system only during fault on and 354 off times. It computes the entire F and then zeros out the part corresponding to 355 differential equations 356 F = [0;g(y)]; 357 */ 358 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx) 359 { 360 PetscErrorCode ierr; 361 Userctx *user=(Userctx*)ctx; 362 PetscScalar *f; 363 PetscInt i; 364 365 PetscFunctionBegin; 366 ierr = ResidualFunction(snes,X,F,user);CHKERRQ(ierr); 367 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 368 for (i=0; i < ngen; i++) { 369 f[9*i] = 0; 370 f[9*i+1] = 0; 371 f[9*i+2] = 0; 372 f[9*i+3] = 0; 373 f[9*i+6] = 0; 374 f[9*i+7] = 0; 375 f[9*i+8] = 0; 376 } 377 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 378 PetscFunctionReturn(0); 379 } 380 381 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user) 382 { 383 PetscErrorCode ierr; 384 PetscInt *d_nnz; 385 PetscInt i,idx=0,start=0; 386 PetscInt ncols; 387 388 PetscFunctionBegin; 389 ierr = PetscMalloc1(user->neqs_pgrid,&d_nnz);CHKERRQ(ierr); 390 for (i=0; i<user->neqs_pgrid; i++) d_nnz[i] = 0; 391 /* Generator subsystem */ 392 for (i=0; i < ngen; i++) { 393 394 d_nnz[idx] += 3; 395 d_nnz[idx+1] += 2; 396 d_nnz[idx+2] += 2; 397 d_nnz[idx+3] += 5; 398 d_nnz[idx+4] += 6; 399 d_nnz[idx+5] += 6; 400 401 d_nnz[user->neqs_gen+2*gbus[i]] += 3; 402 d_nnz[user->neqs_gen+2*gbus[i]+1] += 3; 403 404 d_nnz[idx+6] += 2; 405 d_nnz[idx+7] += 2; 406 d_nnz[idx+8] += 5; 407 408 idx = idx + 9; 409 } 410 411 start = user->neqs_gen; 412 413 for (i=0; i < nbus; i++) { 414 ierr = MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 415 d_nnz[start+2*i] += ncols; 416 d_nnz[start+2*i+1] += ncols; 417 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 418 } 419 420 ierr = MatSeqAIJSetPreallocation(J,0,d_nnz);CHKERRQ(ierr); 421 422 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 423 PetscFunctionReturn(0); 424 } 425 426 /* 427 J = [-df_dx, -df_dy 428 dg_dx, dg_dy] 429 */ 430 PetscErrorCode ResidualJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx) 431 { 432 PetscErrorCode ierr; 433 Userctx *user=(Userctx*)ctx; 434 Vec Xgen,Xnet; 435 PetscScalar *xgen,*xnet; 436 PetscInt i,idx=0; 437 PetscScalar Vr,Vi,Vm,Vm2; 438 PetscScalar Eqp,Edp,delta; /* Generator variables */ 439 PetscScalar Efd; /* Exciter variables */ 440 PetscScalar Id,Iq; /* Generator dq axis currents */ 441 PetscScalar Vd,Vq; 442 PetscScalar val[10]; 443 PetscInt row[2],col[10]; 444 PetscInt net_start=user->neqs_gen; 445 PetscScalar Zdq_inv[4],det; 446 PetscScalar dVd_dVr,dVd_dVi,dVq_dVr,dVq_dVi,dVd_ddelta,dVq_ddelta; 447 PetscScalar dIGr_ddelta,dIGi_ddelta,dIGr_dId,dIGr_dIq,dIGi_dId,dIGi_dIq; 448 PetscScalar dSE_dEfd; 449 PetscScalar dVm_dVd,dVm_dVq,dVm_dVr,dVm_dVi; 450 PetscInt ncols; 451 const PetscInt *cols; 452 const PetscScalar *yvals; 453 PetscInt k; 454 PetscScalar PD,QD,Vm0,*v0,Vm4; 455 PetscScalar dPD_dVr,dPD_dVi,dQD_dVr,dQD_dVi; 456 PetscScalar dIDr_dVr,dIDr_dVi,dIDi_dVr,dIDi_dVi; 457 458 PetscFunctionBegin; 459 ierr = MatZeroEntries(B);CHKERRQ(ierr); 460 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 461 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 462 463 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 464 ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr); 465 466 /* Generator subsystem */ 467 for (i=0; i < ngen; i++) { 468 Eqp = xgen[idx]; 469 Edp = xgen[idx+1]; 470 delta = xgen[idx+2]; 471 Id = xgen[idx+4]; 472 Iq = xgen[idx+5]; 473 Efd = xgen[idx+6]; 474 475 /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */ 476 row[0] = idx; 477 col[0] = idx; col[1] = idx+4; col[2] = idx+6; 478 val[0] = 1/ Td0p[i]; val[1] = (Xd[i] - Xdp[i])/ Td0p[i]; val[2] = -1/Td0p[i]; 479 480 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 481 482 /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */ 483 row[0] = idx + 1; 484 col[0] = idx + 1; col[1] = idx+5; 485 val[0] = 1/Tq0p[i]; val[1] = -(Xq[i] - Xqp[i])/Tq0p[i]; 486 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 487 488 /* fgen[idx+2] = - w + w_s; */ 489 row[0] = idx + 2; 490 col[0] = idx + 2; col[1] = idx + 3; 491 val[0] = 0; val[1] = -1; 492 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 493 494 /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */ 495 row[0] = idx + 3; 496 col[0] = idx; col[1] = idx + 1; col[2] = idx + 3; col[3] = idx + 4; col[4] = idx + 5; 497 val[0] = Iq/M[i]; val[1] = Id/M[i]; val[2] = D[i]/M[i]; val[3] = (Edp + (Xqp[i]-Xdp[i])*Iq)/M[i]; val[4] = (Eqp + (Xqp[i] - Xdp[i])*Id)/M[i]; 498 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 499 500 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 501 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 502 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 503 504 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 505 506 Zdq_inv[0] = Rs[i]/det; 507 Zdq_inv[1] = Xqp[i]/det; 508 Zdq_inv[2] = -Xdp[i]/det; 509 Zdq_inv[3] = Rs[i]/det; 510 511 dVd_dVr = PetscSinScalar(delta); dVd_dVi = -PetscCosScalar(delta); 512 dVq_dVr = PetscCosScalar(delta); dVq_dVi = PetscSinScalar(delta); 513 dVd_ddelta = Vr*PetscCosScalar(delta) + Vi*PetscSinScalar(delta); 514 dVq_ddelta = -Vr*PetscSinScalar(delta) + Vi*PetscCosScalar(delta); 515 516 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */ 517 row[0] = idx+4; 518 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 519 val[0] = -Zdq_inv[1]; val[1] = -Zdq_inv[0]; val[2] = Zdq_inv[0]*dVd_ddelta + Zdq_inv[1]*dVq_ddelta; 520 col[3] = idx + 4; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 521 val[3] = 1; val[4] = Zdq_inv[0]*dVd_dVr + Zdq_inv[1]*dVq_dVr; val[5] = Zdq_inv[0]*dVd_dVi + Zdq_inv[1]*dVq_dVi; 522 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 523 524 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */ 525 row[0] = idx+5; 526 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 527 val[0] = -Zdq_inv[3]; val[1] = -Zdq_inv[2]; val[2] = Zdq_inv[2]*dVd_ddelta + Zdq_inv[3]*dVq_ddelta; 528 col[3] = idx + 5; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 529 val[3] = 1; val[4] = Zdq_inv[2]*dVd_dVr + Zdq_inv[3]*dVq_dVr; val[5] = Zdq_inv[2]*dVd_dVi + Zdq_inv[3]*dVq_dVi; 530 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 531 532 dIGr_ddelta = Id*PetscCosScalar(delta) - Iq*PetscSinScalar(delta); 533 dIGi_ddelta = Id*PetscSinScalar(delta) + Iq*PetscCosScalar(delta); 534 dIGr_dId = PetscSinScalar(delta); dIGr_dIq = PetscCosScalar(delta); 535 dIGi_dId = -PetscCosScalar(delta); dIGi_dIq = PetscSinScalar(delta); 536 537 /* fnet[2*gbus[i]] -= IGi; */ 538 row[0] = net_start + 2*gbus[i]; 539 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 540 val[0] = -dIGi_ddelta; val[1] = -dIGi_dId; val[2] = -dIGi_dIq; 541 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 542 543 /* fnet[2*gbus[i]+1] -= IGr; */ 544 row[0] = net_start + 2*gbus[i]+1; 545 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 546 val[0] = -dIGr_ddelta; val[1] = -dIGr_dId; val[2] = -dIGr_dIq; 547 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 548 549 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 550 551 /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */ 552 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */ 553 dSE_dEfd = k1[i]*k2[i]*PetscExpScalar(k2[i]*Efd); 554 555 row[0] = idx + 6; 556 col[0] = idx + 6; col[1] = idx + 8; 557 val[0] = (KE[i] + dSE_dEfd)/TE[i]; val[1] = -1/TE[i]; 558 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 559 560 /* Exciter differential equations */ 561 562 /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */ 563 row[0] = idx + 7; 564 col[0] = idx + 6; col[1] = idx + 7; 565 val[0] = (-KF[i]/TF[i])/TF[i]; val[1] = 1/TF[i]; 566 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 567 568 /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */ 569 /* Vm = (Vd^2 + Vq^2)^0.5; */ 570 dVm_dVd = Vd/Vm; dVm_dVq = Vq/Vm; 571 dVm_dVr = dVm_dVd*dVd_dVr + dVm_dVq*dVq_dVr; 572 dVm_dVi = dVm_dVd*dVd_dVi + dVm_dVq*dVq_dVi; 573 row[0] = idx + 8; 574 col[0] = idx + 6; col[1] = idx + 7; col[2] = idx + 8; 575 val[0] = (KA[i]*KF[i]/TF[i])/TA[i]; val[1] = -KA[i]/TA[i]; val[2] = 1/TA[i]; 576 col[3] = net_start + 2*gbus[i]; col[4] = net_start + 2*gbus[i]+1; 577 val[3] = KA[i]*dVm_dVr/TA[i]; val[4] = KA[i]*dVm_dVi/TA[i]; 578 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 579 idx = idx + 9; 580 } 581 582 for (i=0; i<nbus; i++) { 583 ierr = MatGetRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 584 row[0] = net_start + 2*i; 585 for (k=0; k<ncols; k++) { 586 col[k] = net_start + cols[k]; 587 val[k] = yvals[k]; 588 } 589 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 590 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 591 592 ierr = MatGetRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 593 row[0] = net_start + 2*i+1; 594 for (k=0; k<ncols; k++) { 595 col[k] = net_start + cols[k]; 596 val[k] = yvals[k]; 597 } 598 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 599 ierr = MatRestoreRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 600 } 601 602 ierr = MatAssemblyBegin(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 603 ierr = MatAssemblyEnd(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 604 605 606 ierr = VecGetArray(user->V0,&v0);CHKERRQ(ierr); 607 for (i=0; i < nload; i++) { 608 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 609 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 610 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; Vm4 = Vm2*Vm2; 611 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 612 PD = QD = 0.0; 613 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0; 614 for (k=0; k < ld_nsegsp[i]; k++) { 615 PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 616 dPD_dVr += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vr*PetscPowScalar(Vm,(ld_betap[k]-2)); 617 dPD_dVi += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vi*PetscPowScalar(Vm,(ld_betap[k]-2)); 618 } 619 for (k=0; k < ld_nsegsq[i]; k++) { 620 QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 621 dQD_dVr += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vr*PetscPowScalar(Vm,(ld_betaq[k]-2)); 622 dQD_dVi += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vi*PetscPowScalar(Vm,(ld_betaq[k]-2)); 623 } 624 625 /* IDr = (PD*Vr + QD*Vi)/Vm2; */ 626 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */ 627 628 dIDr_dVr = (dPD_dVr*Vr + dQD_dVr*Vi + PD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vr)/Vm4; 629 dIDr_dVi = (dPD_dVi*Vr + dQD_dVi*Vi + QD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vi)/Vm4; 630 631 dIDi_dVr = (-dQD_dVr*Vr + dPD_dVr*Vi - QD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vr)/Vm4; 632 dIDi_dVi = (-dQD_dVi*Vr + dPD_dVi*Vi + PD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vi)/Vm4; 633 634 635 /* fnet[2*lbus[i]] += IDi; */ 636 row[0] = net_start + 2*lbus[i]; 637 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 638 val[0] = dIDi_dVr; val[1] = dIDi_dVi; 639 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 640 /* fnet[2*lbus[i]+1] += IDr; */ 641 row[0] = net_start + 2*lbus[i]+1; 642 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 643 val[0] = dIDr_dVr; val[1] = dIDr_dVi; 644 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 645 } 646 ierr = VecRestoreArray(user->V0,&v0);CHKERRQ(ierr); 647 648 ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr); 649 ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr); 650 651 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 652 653 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 654 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 655 PetscFunctionReturn(0); 656 } 657 658 /* 659 J = [I, 0 660 dg_dx, dg_dy] 661 */ 662 PetscErrorCode AlgJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx) 663 { 664 PetscErrorCode ierr; 665 Userctx *user=(Userctx*)ctx; 666 667 PetscFunctionBegin; 668 ierr = ResidualJacobian(snes,X,A,B,ctx);CHKERRQ(ierr); 669 ierr = MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);CHKERRQ(ierr); 670 ierr = MatZeroRowsIS(A,user->is_diff,1.0,NULL,NULL);CHKERRQ(ierr); 671 PetscFunctionReturn(0); 672 } 673 674 /* 675 J = [a*I-df_dx, -df_dy 676 dg_dx, dg_dy] 677 */ 678 679 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx *user) 680 { 681 PetscErrorCode ierr; 682 SNES snes; 683 PetscScalar atmp = (PetscScalar) a; 684 PetscInt i,row; 685 686 PetscFunctionBegin; 687 user->t = t; 688 689 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 690 ierr = ResidualJacobian(snes,X,A,B,user);CHKERRQ(ierr); 691 for (i=0;i < ngen;i++) { 692 row = 9*i; 693 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 694 row = 9*i+1; 695 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 696 row = 9*i+2; 697 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 698 row = 9*i+3; 699 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 700 row = 9*i+6; 701 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 702 row = 9*i+7; 703 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 704 row = 9*i+8; 705 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 706 } 707 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 708 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 709 PetscFunctionReturn(0); 710 } 711 712 int main(int argc,char **argv) 713 { 714 TS ts; 715 SNES snes_alg; 716 PetscErrorCode ierr; 717 PetscMPIInt size; 718 Userctx user; 719 PetscViewer Xview,Ybusview; 720 Vec X; 721 Mat J; 722 PetscInt i; 723 /* sensitivity context */ 724 PetscScalar *y_ptr; 725 Vec lambda[1]; 726 PetscInt *idx2; 727 Vec Xdot; 728 Vec F_alg; 729 PetscInt row_loc,col_loc; 730 PetscScalar val; 731 732 ierr = PetscInitialize(&argc,&argv,"petscoptions",help);if (ierr) return ierr; 733 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 734 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 735 736 user.neqs_gen = 9*ngen; /* # eqs. for generator subsystem */ 737 user.neqs_net = 2*nbus; /* # eqs. for network subsystem */ 738 user.neqs_pgrid = user.neqs_gen + user.neqs_net; 739 740 /* Create indices for differential and algebraic equations */ 741 ierr = PetscMalloc1(7*ngen,&idx2);CHKERRQ(ierr); 742 for (i=0; i<ngen; i++) { 743 idx2[7*i] = 9*i; idx2[7*i+1] = 9*i+1; idx2[7*i+2] = 9*i+2; idx2[7*i+3] = 9*i+3; 744 idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8; 745 } 746 ierr = ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff);CHKERRQ(ierr); 747 ierr = ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg);CHKERRQ(ierr); 748 ierr = PetscFree(idx2);CHKERRQ(ierr); 749 750 /* Read initial voltage vector and Ybus */ 751 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"X.bin",FILE_MODE_READ,&Xview);CHKERRQ(ierr); 752 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ybus.bin",FILE_MODE_READ,&Ybusview);CHKERRQ(ierr); 753 754 ierr = VecCreate(PETSC_COMM_WORLD,&user.V0);CHKERRQ(ierr); 755 ierr = VecSetSizes(user.V0,PETSC_DECIDE,user.neqs_net);CHKERRQ(ierr); 756 ierr = VecLoad(user.V0,Xview);CHKERRQ(ierr); 757 758 ierr = MatCreate(PETSC_COMM_WORLD,&user.Ybus);CHKERRQ(ierr); 759 ierr = MatSetSizes(user.Ybus,PETSC_DECIDE,PETSC_DECIDE,user.neqs_net,user.neqs_net);CHKERRQ(ierr); 760 ierr = MatSetType(user.Ybus,MATBAIJ);CHKERRQ(ierr); 761 /* ierr = MatSetBlockSize(user.Ybus,2);CHKERRQ(ierr); */ 762 ierr = MatLoad(user.Ybus,Ybusview);CHKERRQ(ierr); 763 764 /* Set run time options */ 765 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");CHKERRQ(ierr); 766 { 767 user.tfaulton = 1.0; 768 user.tfaultoff = 1.2; 769 user.Rfault = 0.0001; 770 user.faultbus = 8; 771 ierr = PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL);CHKERRQ(ierr); 772 ierr = PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL);CHKERRQ(ierr); 773 ierr = PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL);CHKERRQ(ierr); 774 user.t0 = 0.0; 775 user.tmax = 5.0; 776 ierr = PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL);CHKERRQ(ierr); 777 ierr = PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL);CHKERRQ(ierr); 778 } 779 ierr = PetscOptionsEnd();CHKERRQ(ierr); 780 781 ierr = PetscViewerDestroy(&Xview);CHKERRQ(ierr); 782 ierr = PetscViewerDestroy(&Ybusview);CHKERRQ(ierr); 783 784 /* Create DMs for generator and network subsystems */ 785 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); 786 ierr = DMSetOptionsPrefix(user.dmgen,"dmgen_");CHKERRQ(ierr); 787 ierr = DMSetFromOptions(user.dmgen);CHKERRQ(ierr); 788 ierr = DMSetUp(user.dmgen);CHKERRQ(ierr); 789 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); 790 ierr = DMSetOptionsPrefix(user.dmnet,"dmnet_");CHKERRQ(ierr); 791 ierr = DMSetFromOptions(user.dmnet);CHKERRQ(ierr); 792 ierr = DMSetUp(user.dmnet);CHKERRQ(ierr); 793 /* Create a composite DM packer and add the two DMs */ 794 ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid);CHKERRQ(ierr); 795 ierr = DMSetOptionsPrefix(user.dmpgrid,"pgrid_");CHKERRQ(ierr); 796 ierr = DMCompositeAddDM(user.dmpgrid,user.dmgen);CHKERRQ(ierr); 797 ierr = DMCompositeAddDM(user.dmpgrid,user.dmnet);CHKERRQ(ierr); 798 799 ierr = DMCreateGlobalVector(user.dmpgrid,&X);CHKERRQ(ierr); 800 801 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 802 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,user.neqs_pgrid,user.neqs_pgrid);CHKERRQ(ierr); 803 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 804 ierr = PreallocateJacobian(J,&user);CHKERRQ(ierr); 805 806 807 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 808 Create timestepping solver context 809 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 810 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 811 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 812 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 813 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&user);CHKERRQ(ierr); 814 ierr = TSSetIJacobian(ts,J,J,(TSIJacobian)IJacobian,&user);CHKERRQ(ierr); 815 ierr = TSSetApplicationContext(ts,&user);CHKERRQ(ierr); 816 817 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 818 Set initial conditions 819 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 820 ierr = SetInitialGuess(X,&user);CHKERRQ(ierr); 821 /* Just to set up the Jacobian structure */ 822 ierr = VecDuplicate(X,&Xdot);CHKERRQ(ierr); 823 ierr = IJacobian(ts,0.0,X,Xdot,0.0,J,J,&user);CHKERRQ(ierr); 824 ierr = VecDestroy(&Xdot);CHKERRQ(ierr); 825 826 /* 827 Save trajectory of solution so that TSAdjointSolve() may be used 828 */ 829 ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr); 830 831 ierr = TSSetMaxTime(ts,user.tmax);CHKERRQ(ierr); 832 ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr); 833 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 834 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 835 836 user.alg_flg = PETSC_FALSE; 837 /* Prefault period */ 838 ierr = TSSolve(ts,X);CHKERRQ(ierr); 839 840 /* Create the nonlinear solver for solving the algebraic system */ 841 /* Note that although the algebraic system needs to be solved only for 842 Idq and V, we reuse the entire system including xgen. The xgen 843 variables are held constant by setting their residuals to 0 and 844 putting a 1 on the Jacobian diagonal for xgen rows 845 */ 846 ierr = VecDuplicate(X,&F_alg);CHKERRQ(ierr); 847 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_alg);CHKERRQ(ierr); 848 ierr = SNESSetFunction(snes_alg,F_alg,AlgFunction,&user);CHKERRQ(ierr); 849 ierr = MatZeroEntries(J);CHKERRQ(ierr); 850 ierr = SNESSetJacobian(snes_alg,J,J,AlgJacobian,&user);CHKERRQ(ierr); 851 ierr = SNESSetOptionsPrefix(snes_alg,"alg_");CHKERRQ(ierr); 852 ierr = SNESSetFromOptions(snes_alg);CHKERRQ(ierr); 853 854 /* Apply disturbance - resistive fault at user.faultbus */ 855 /* This is done by adding shunt conductance to the diagonal location 856 in the Ybus matrix */ 857 row_loc = 2*user.faultbus; col_loc = 2*user.faultbus+1; /* Location for G */ 858 val = 1/user.Rfault; 859 ierr = MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 860 row_loc = 2*user.faultbus+1; col_loc = 2*user.faultbus; /* Location for G */ 861 val = 1/user.Rfault; 862 ierr = MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 863 864 ierr = MatAssemblyBegin(user.Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 865 ierr = MatAssemblyEnd(user.Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 866 867 user.alg_flg = PETSC_TRUE; 868 /* Solve the algebraic equations */ 869 ierr = SNESSolve(snes_alg,NULL,X);CHKERRQ(ierr); 870 871 872 /* Disturbance period */ 873 user.alg_flg = PETSC_FALSE; 874 ierr = TSSetTime(ts,user.tfaulton);CHKERRQ(ierr); 875 ierr = TSSetMaxTime(ts,user.tfaultoff);CHKERRQ(ierr); 876 ierr = TSSolve(ts,X);CHKERRQ(ierr); 877 878 /* Remove the fault */ 879 row_loc = 2*user.faultbus; col_loc = 2*user.faultbus+1; 880 val = -1/user.Rfault; 881 ierr = MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 882 row_loc = 2*user.faultbus+1; col_loc = 2*user.faultbus; 883 val = -1/user.Rfault; 884 ierr = MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 885 886 ierr = MatAssemblyBegin(user.Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 887 ierr = MatAssemblyEnd(user.Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 888 889 ierr = MatZeroEntries(J);CHKERRQ(ierr); 890 891 user.alg_flg = PETSC_TRUE; 892 893 /* Solve the algebraic equations */ 894 ierr = SNESSolve(snes_alg,NULL,X);CHKERRQ(ierr); 895 896 /* Post-disturbance period */ 897 user.alg_flg = PETSC_TRUE; 898 ierr = TSSetTime(ts,user.tfaultoff);CHKERRQ(ierr); 899 ierr = TSSetMaxTime(ts,user.tmax);CHKERRQ(ierr); 900 ierr = TSSolve(ts,X);CHKERRQ(ierr); 901 902 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 903 Adjoint model starts here 904 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 905 ierr = TSSetPostStep(ts,NULL);CHKERRQ(ierr); 906 ierr = MatCreateVecs(J,&lambda[0],NULL);CHKERRQ(ierr); 907 /* Set initial conditions for the adjoint integration */ 908 ierr = VecZeroEntries(lambda[0]);CHKERRQ(ierr); 909 ierr = VecGetArray(lambda[0],&y_ptr);CHKERRQ(ierr); 910 y_ptr[0] = 1.0; 911 ierr = VecRestoreArray(lambda[0],&y_ptr);CHKERRQ(ierr); 912 ierr = TSSetCostGradients(ts,1,lambda,NULL);CHKERRQ(ierr); 913 914 ierr = TSAdjointSolve(ts);CHKERRQ(ierr); 915 916 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: \n");CHKERRQ(ierr); 917 ierr = VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 918 ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr); 919 920 ierr = SNESDestroy(&snes_alg);CHKERRQ(ierr); 921 ierr = VecDestroy(&F_alg);CHKERRQ(ierr); 922 ierr = MatDestroy(&J);CHKERRQ(ierr); 923 ierr = MatDestroy(&user.Ybus);CHKERRQ(ierr); 924 ierr = VecDestroy(&X);CHKERRQ(ierr); 925 ierr = VecDestroy(&user.V0);CHKERRQ(ierr); 926 ierr = DMDestroy(&user.dmgen);CHKERRQ(ierr); 927 ierr = DMDestroy(&user.dmnet);CHKERRQ(ierr); 928 ierr = DMDestroy(&user.dmpgrid);CHKERRQ(ierr); 929 ierr = ISDestroy(&user.is_diff);CHKERRQ(ierr); 930 ierr = ISDestroy(&user.is_alg);CHKERRQ(ierr); 931 ierr = TSDestroy(&ts);CHKERRQ(ierr); 932 ierr = PetscFinalize(); 933 return ierr; 934 } 935 936 /*TEST 937 938 build: 939 requires: double !complex !define(PETSC_USE_64BIT_INDICES) 940 941 test: 942 args: -viewer_binary_skip_info 943 localrunfiles: petscoptions X.bin Ybus.bin 944 945 TEST*/ 946