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