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 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 131 132 /* Network subsystem initialization */ 133 PetscCall(VecCopy(user->V0,Xnet)); 134 135 /* Generator subsystem initialization */ 136 PetscCall(VecGetArray(Xgen,&xgen)); 137 PetscCall(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 PetscCall(VecRestoreArray(Xgen,&xgen)); 190 PetscCall(VecRestoreArray(Xnet,&xnet)); 191 192 /* PetscCall(VecView(Xgen,0)); */ 193 PetscCall(DMCompositeGather(user->dmpgrid,INSERT_VALUES,X,Xgen,Xnet)); 194 PetscCall(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 PetscCall(VecZeroEntries(F)); 216 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 217 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid,&Fgen,&Fnet)); 218 PetscCall(DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet)); 219 PetscCall(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 PetscCall(MatMult(user->Ybus,Xnet,Fnet)); 230 231 PetscCall(VecGetArray(Xgen,&xgen)); 232 PetscCall(VecGetArray(Xnet,&xnet)); 233 PetscCall(VecGetArray(Fgen,&fgen)); 234 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(VecRestoreArray(user->V0,&v0)); 306 307 PetscCall(VecRestoreArray(Xgen,&xgen)); 308 PetscCall(VecRestoreArray(Xnet,&xnet)); 309 PetscCall(VecRestoreArray(Fgen,&fgen)); 310 PetscCall(VecRestoreArray(Fnet,&fnet)); 311 312 PetscCall(DMCompositeGather(user->dmpgrid,INSERT_VALUES,F,Fgen,Fnet)); 313 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 314 PetscCall(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 PetscCall(TSGetSNES(ts,&snes)); 332 PetscCall(ResidualFunction(snes,X,F,user)); 333 PetscCall(VecGetArray(F,&f)); 334 PetscCall(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 PetscCall(VecRestoreArray(F,&f)); 345 PetscCall(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 PetscCall(ResidualFunction(snes,X,F,user)); 362 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL)); 409 d_nnz[start+2*i] += ncols; 410 d_nnz[start+2*i+1] += ncols; 411 PetscCall(MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL)); 412 } 413 414 PetscCall(MatSeqAIJSetPreallocation(J,0,d_nnz)); 415 416 PetscCall(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 PetscCall(MatZeroEntries(B)); 453 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 454 PetscCall(DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet)); 455 456 PetscCall(VecGetArray(Xgen,&xgen)); 457 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(J,1,row,5,col,val,INSERT_VALUES)); 572 idx = idx + 9; 573 } 574 575 for (i=0; i<nbus; i++) { 576 PetscCall(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 PetscCall(MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES)); 583 PetscCall(MatRestoreRow(user->Ybus,2*i,&ncols,&cols,&yvals)); 584 585 PetscCall(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 PetscCall(MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES)); 592 PetscCall(MatRestoreRow(user->Ybus,2*i+1,&ncols,&cols,&yvals)); 593 } 594 595 PetscCall(MatAssemblyBegin(J,MAT_FLUSH_ASSEMBLY)); 596 PetscCall(MatAssemblyEnd(J,MAT_FLUSH_ASSEMBLY)); 597 598 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(J,1,row,2,col,val,ADD_VALUES)); 636 } 637 PetscCall(VecRestoreArray(user->V0,&v0)); 638 639 PetscCall(VecRestoreArray(Xgen,&xgen)); 640 PetscCall(VecRestoreArray(Xnet,&xnet)); 641 642 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet)); 643 644 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 645 PetscCall(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 PetscCall(ResidualJacobian(snes,X,A,B,ctx)); 659 PetscCall(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE)); 660 PetscCall(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 PetscCall(TSGetSNES(ts,&snes)); 679 PetscCall(ResidualJacobian(snes,X,A,B,user)); 680 for (i=0;i < ngen;i++) { 681 row = 9*i; 682 PetscCall(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 683 row = 9*i+1; 684 PetscCall(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 685 row = 9*i+2; 686 PetscCall(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 687 row = 9*i+3; 688 PetscCall(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 689 row = 9*i+6; 690 PetscCall(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 691 row = 9*i+7; 692 PetscCall(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 693 row = 9*i+8; 694 PetscCall(MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES)); 695 } 696 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 697 PetscCall(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 PetscMPIInt size; 706 Userctx user; 707 PetscViewer Xview,Ybusview; 708 Vec X; 709 Mat J; 710 PetscInt i; 711 /* sensitivity context */ 712 PetscScalar *y_ptr; 713 Vec lambda[1]; 714 PetscInt *idx2; 715 Vec Xdot; 716 Vec F_alg; 717 PetscInt row_loc,col_loc; 718 PetscScalar val; 719 720 PetscCall(PetscInitialize(&argc,&argv,"petscoptions",help)); 721 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 722 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 723 724 user.neqs_gen = 9*ngen; /* # eqs. for generator subsystem */ 725 user.neqs_net = 2*nbus; /* # eqs. for network subsystem */ 726 user.neqs_pgrid = user.neqs_gen + user.neqs_net; 727 728 /* Create indices for differential and algebraic equations */ 729 PetscCall(PetscMalloc1(7*ngen,&idx2)); 730 for (i=0; i<ngen; i++) { 731 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; 732 idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8; 733 } 734 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff)); 735 PetscCall(ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg)); 736 PetscCall(PetscFree(idx2)); 737 738 /* Read initial voltage vector and Ybus */ 739 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"X.bin",FILE_MODE_READ,&Xview)); 740 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ybus.bin",FILE_MODE_READ,&Ybusview)); 741 742 PetscCall(VecCreate(PETSC_COMM_WORLD,&user.V0)); 743 PetscCall(VecSetSizes(user.V0,PETSC_DECIDE,user.neqs_net)); 744 PetscCall(VecLoad(user.V0,Xview)); 745 746 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.Ybus)); 747 PetscCall(MatSetSizes(user.Ybus,PETSC_DECIDE,PETSC_DECIDE,user.neqs_net,user.neqs_net)); 748 PetscCall(MatSetType(user.Ybus,MATBAIJ)); 749 /* PetscCall(MatSetBlockSize(user.Ybus,2)); */ 750 PetscCall(MatLoad(user.Ybus,Ybusview)); 751 752 /* Set run time options */ 753 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options",""); 754 { 755 user.tfaulton = 1.0; 756 user.tfaultoff = 1.2; 757 user.Rfault = 0.0001; 758 user.faultbus = 8; 759 PetscCall(PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL)); 760 PetscCall(PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL)); 761 PetscCall(PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL)); 762 user.t0 = 0.0; 763 user.tmax = 5.0; 764 PetscCall(PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL)); 765 PetscCall(PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL)); 766 } 767 PetscOptionsEnd(); 768 769 PetscCall(PetscViewerDestroy(&Xview)); 770 PetscCall(PetscViewerDestroy(&Ybusview)); 771 772 /* Create DMs for generator and network subsystems */ 773 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen)); 774 PetscCall(DMSetOptionsPrefix(user.dmgen,"dmgen_")); 775 PetscCall(DMSetFromOptions(user.dmgen)); 776 PetscCall(DMSetUp(user.dmgen)); 777 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet)); 778 PetscCall(DMSetOptionsPrefix(user.dmnet,"dmnet_")); 779 PetscCall(DMSetFromOptions(user.dmnet)); 780 PetscCall(DMSetUp(user.dmnet)); 781 /* Create a composite DM packer and add the two DMs */ 782 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid)); 783 PetscCall(DMSetOptionsPrefix(user.dmpgrid,"pgrid_")); 784 PetscCall(DMCompositeAddDM(user.dmpgrid,user.dmgen)); 785 PetscCall(DMCompositeAddDM(user.dmpgrid,user.dmnet)); 786 787 PetscCall(DMCreateGlobalVector(user.dmpgrid,&X)); 788 789 PetscCall(MatCreate(PETSC_COMM_WORLD,&J)); 790 PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,user.neqs_pgrid,user.neqs_pgrid)); 791 PetscCall(MatSetFromOptions(J)); 792 PetscCall(PreallocateJacobian(J,&user)); 793 794 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 795 Create timestepping solver context 796 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 797 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 798 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 799 PetscCall(TSSetType(ts,TSCN)); 800 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&user)); 801 PetscCall(TSSetIJacobian(ts,J,J,(TSIJacobian)IJacobian,&user)); 802 PetscCall(TSSetApplicationContext(ts,&user)); 803 804 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 805 Set initial conditions 806 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 807 PetscCall(SetInitialGuess(X,&user)); 808 /* Just to set up the Jacobian structure */ 809 PetscCall(VecDuplicate(X,&Xdot)); 810 PetscCall(IJacobian(ts,0.0,X,Xdot,0.0,J,J,&user)); 811 PetscCall(VecDestroy(&Xdot)); 812 813 /* 814 Save trajectory of solution so that TSAdjointSolve() may be used 815 */ 816 PetscCall(TSSetSaveTrajectory(ts)); 817 818 PetscCall(TSSetMaxTime(ts,user.tmax)); 819 PetscCall(TSSetTimeStep(ts,0.01)); 820 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 821 PetscCall(TSSetFromOptions(ts)); 822 823 user.alg_flg = PETSC_FALSE; 824 /* Prefault period */ 825 PetscCall(TSSolve(ts,X)); 826 827 /* Create the nonlinear solver for solving the algebraic system */ 828 /* Note that although the algebraic system needs to be solved only for 829 Idq and V, we reuse the entire system including xgen. The xgen 830 variables are held constant by setting their residuals to 0 and 831 putting a 1 on the Jacobian diagonal for xgen rows 832 */ 833 PetscCall(VecDuplicate(X,&F_alg)); 834 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes_alg)); 835 PetscCall(SNESSetFunction(snes_alg,F_alg,AlgFunction,&user)); 836 PetscCall(MatZeroEntries(J)); 837 PetscCall(SNESSetJacobian(snes_alg,J,J,AlgJacobian,&user)); 838 PetscCall(SNESSetOptionsPrefix(snes_alg,"alg_")); 839 PetscCall(SNESSetFromOptions(snes_alg)); 840 841 /* Apply disturbance - resistive fault at user.faultbus */ 842 /* This is done by adding shunt conductance to the diagonal location 843 in the Ybus matrix */ 844 row_loc = 2*user.faultbus; col_loc = 2*user.faultbus+1; /* Location for G */ 845 val = 1/user.Rfault; 846 PetscCall(MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 847 row_loc = 2*user.faultbus+1; col_loc = 2*user.faultbus; /* Location for G */ 848 val = 1/user.Rfault; 849 PetscCall(MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 850 851 PetscCall(MatAssemblyBegin(user.Ybus,MAT_FINAL_ASSEMBLY)); 852 PetscCall(MatAssemblyEnd(user.Ybus,MAT_FINAL_ASSEMBLY)); 853 854 user.alg_flg = PETSC_TRUE; 855 /* Solve the algebraic equations */ 856 PetscCall(SNESSolve(snes_alg,NULL,X)); 857 858 /* Disturbance period */ 859 user.alg_flg = PETSC_FALSE; 860 PetscCall(TSSetTime(ts,user.tfaulton)); 861 PetscCall(TSSetMaxTime(ts,user.tfaultoff)); 862 PetscCall(TSSolve(ts,X)); 863 864 /* Remove the fault */ 865 row_loc = 2*user.faultbus; col_loc = 2*user.faultbus+1; 866 val = -1/user.Rfault; 867 PetscCall(MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 868 row_loc = 2*user.faultbus+1; col_loc = 2*user.faultbus; 869 val = -1/user.Rfault; 870 PetscCall(MatSetValues(user.Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES)); 871 872 PetscCall(MatAssemblyBegin(user.Ybus,MAT_FINAL_ASSEMBLY)); 873 PetscCall(MatAssemblyEnd(user.Ybus,MAT_FINAL_ASSEMBLY)); 874 875 PetscCall(MatZeroEntries(J)); 876 877 user.alg_flg = PETSC_TRUE; 878 879 /* Solve the algebraic equations */ 880 PetscCall(SNESSolve(snes_alg,NULL,X)); 881 882 /* Post-disturbance period */ 883 user.alg_flg = PETSC_TRUE; 884 PetscCall(TSSetTime(ts,user.tfaultoff)); 885 PetscCall(TSSetMaxTime(ts,user.tmax)); 886 PetscCall(TSSolve(ts,X)); 887 888 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 889 Adjoint model starts here 890 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 891 PetscCall(TSSetPostStep(ts,NULL)); 892 PetscCall(MatCreateVecs(J,&lambda[0],NULL)); 893 /* Set initial conditions for the adjoint integration */ 894 PetscCall(VecZeroEntries(lambda[0])); 895 PetscCall(VecGetArray(lambda[0],&y_ptr)); 896 y_ptr[0] = 1.0; 897 PetscCall(VecRestoreArray(lambda[0],&y_ptr)); 898 PetscCall(TSSetCostGradients(ts,1,lambda,NULL)); 899 900 PetscCall(TSAdjointSolve(ts)); 901 902 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: \n")); 903 PetscCall(VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD)); 904 PetscCall(VecDestroy(&lambda[0])); 905 906 PetscCall(SNESDestroy(&snes_alg)); 907 PetscCall(VecDestroy(&F_alg)); 908 PetscCall(MatDestroy(&J)); 909 PetscCall(MatDestroy(&user.Ybus)); 910 PetscCall(VecDestroy(&X)); 911 PetscCall(VecDestroy(&user.V0)); 912 PetscCall(DMDestroy(&user.dmgen)); 913 PetscCall(DMDestroy(&user.dmnet)); 914 PetscCall(DMDestroy(&user.dmpgrid)); 915 PetscCall(ISDestroy(&user.is_diff)); 916 PetscCall(ISDestroy(&user.is_alg)); 917 PetscCall(TSDestroy(&ts)); 918 PetscCall(PetscFinalize()); 919 return 0; 920 } 921 922 /*TEST 923 924 build: 925 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) 926 927 test: 928 args: -viewer_binary_skip_info 929 localrunfiles: petscoptions X.bin Ybus.bin 930 931 test: 932 suffix: 2 933 args: -viewer_binary_skip_info -ts_type beuler 934 localrunfiles: petscoptions X.bin Ybus.bin 935 936 TEST*/ 937