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