1 2 static char help[] = "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 11 12 \dot{x} = f(x,y,t) 13 0 = g(x,y,t) 14 15 where the generators are described by differential equations, while the algebraic 16 constraints define the network equations. 17 18 The generators are modeled with a 4th order differential equation describing the electrical 19 and mechanical dynamics. Each generator also has an exciter system modeled by 3rd order 20 diff. eqns. describing the exciter, voltage regulator, and the feedback stabilizer 21 mechanism. 22 23 The network equations are described by nodal current balance equations. 24 I(x,y) - Y*V = 0 25 26 where: 27 I(x,y) is the current injected from generators and loads. 28 Y is the admittance matrix, and 29 V is the voltage vector 30 */ 31 32 /* 33 Include "petscts.h" so that we can use TS solvers. Note that this 34 file automatically includes: 35 petscsys.h - base PETSc routines petscvec.h - vectors 36 petscmat.h - matrices 37 petscis.h - index sets petscksp.h - Krylov subspace methods 38 petscviewer.h - viewers petscpc.h - preconditioners 39 petscksp.h - linear solvers 40 */ 41 42 #include <petscts.h> 43 #include <petscdm.h> 44 #include <petscdmda.h> 45 #include <petscdmcomposite.h> 46 47 #define freq 60 48 #define w_s (2*PETSC_PI*freq) 49 50 /* Sizes and indices */ 51 const PetscInt nbus = 9; /* Number of network buses */ 52 const PetscInt ngen = 3; /* Number of generators */ 53 const PetscInt nload = 3; /* Number of loads */ 54 const PetscInt gbus[3] = {0,1,2}; /* Buses at which generators are incident */ 55 const PetscInt lbus[3] = {4,5,7}; /* Buses at which loads are incident */ 56 57 /* Generator real and reactive powers (found via loadflow) */ 58 const PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000}; 59 const PetscScalar QG[3] = {0.270702180178785,0.066120127797275,-0.108402221791588}; 60 /* Generator constants */ 61 const PetscScalar H[3] = {23.64,6.4,3.01}; /* Inertia constant */ 62 const PetscScalar Rs[3] = {0.0,0.0,0.0}; /* Stator Resistance */ 63 const PetscScalar Xd[3] = {0.146,0.8958,1.3125}; /* d-axis reactance */ 64 const PetscScalar Xdp[3] = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */ 65 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 */ 66 const PetscScalar Xqp[3] = {0.0969,0.1969,0.25}; /* q-axis transient reactance */ 67 const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */ 68 const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */ 69 PetscScalar M[3]; /* M = 2*H/w_s */ 70 PetscScalar D[3]; /* D = 0.1*M */ 71 72 PetscScalar TM[3]; /* Mechanical Torque */ 73 /* Exciter system constants */ 74 const PetscScalar KA[3] = {20.0,20.0,20.0}; /* Voltage regulartor gain constant */ 75 const PetscScalar TA[3] = {0.2,0.2,0.2}; /* Voltage regulator time constant */ 76 const PetscScalar KE[3] = {1.0,1.0,1.0}; /* Exciter gain constant */ 77 const PetscScalar TE[3] = {0.314,0.314,0.314}; /* Exciter time constant */ 78 const PetscScalar KF[3] = {0.063,0.063,0.063}; /* Feedback stabilizer gain constant */ 79 const PetscScalar TF[3] = {0.35,0.35,0.35}; /* Feedback stabilizer time constant */ 80 const PetscScalar k1[3] = {0.0039,0.0039,0.0039}; 81 const PetscScalar k2[3] = {1.555,1.555,1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */ 82 const PetscScalar VRMIN[3] = {-4.0,-4.0,-4.0}; 83 const PetscScalar VRMAX[3] = {7.0,7.0,7.0}; 84 PetscInt VRatmin[3]; 85 PetscInt VRatmax[3]; 86 87 PetscScalar Vref[3]; 88 /* Load constants 89 We use a composite load model that describes the load and reactive powers at each time instant as follows 90 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i 91 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i 92 where 93 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads 94 ld_alphap,ld_alphap - Percentage contribution (weights) or loads 95 P_D0 - Real power load 96 Q_D0 - Reactive power load 97 V_m(t) - Voltage magnitude at time t 98 V_m0 - Voltage magnitude at t = 0 99 ld_betap, ld_betaq - exponents describing the load model for real and reactive part 100 101 Note: All loads have the same characteristic currently. 102 */ 103 const PetscScalar PD0[3] = {1.25,0.9,1.0}; 104 const PetscScalar QD0[3] = {0.5,0.3,0.35}; 105 const PetscInt ld_nsegsp[3] = {3,3,3}; 106 const PetscScalar ld_alphap[3] = {1.0,0.0,0.0}; 107 const PetscScalar ld_betap[3] = {2.0,1.0,0.0}; 108 const PetscInt ld_nsegsq[3] = {3,3,3}; 109 const PetscScalar ld_alphaq[3] = {1.0,0.0,0.0}; 110 const PetscScalar ld_betaq[3] = {2.0,1.0,0.0}; 111 112 typedef struct { 113 DM dmgen, dmnet; /* DMs to manage generator and network subsystem */ 114 DM dmpgrid; /* Composite DM to manage the entire power grid */ 115 Mat Ybus; /* Network admittance matrix */ 116 Vec V0; /* Initial voltage vector (Power flow solution) */ 117 PetscReal tfaulton,tfaultoff; /* Fault on and off times */ 118 PetscInt faultbus; /* Fault bus */ 119 PetscScalar Rfault; 120 PetscReal t0,tmax; 121 PetscInt neqs_gen,neqs_net,neqs_pgrid; 122 Mat Sol; /* Matrix to save solution at each time step */ 123 PetscInt stepnum; 124 PetscReal t; 125 SNES snes_alg; 126 IS is_diff; /* indices for differential equations */ 127 IS is_alg; /* indices for algebraic equations */ 128 PetscBool setisdiff; /* TS computes truncation error based only on the differential variables */ 129 PetscBool semiexplicit; /* If the flag is set then a semi-explicit method is used using TSRK */ 130 } Userctx; 131 132 /* 133 The first two events are for fault on and off, respectively. The following events are 134 to check the min/max limits on the state variable VR. A non windup limiter is used for 135 the VR limits. 136 */ 137 PetscErrorCode EventFunction(TS ts,PetscReal t,Vec X,PetscScalar *fvalue,void *ctx) 138 { 139 Userctx *user=(Userctx*)ctx; 140 Vec Xgen,Xnet; 141 PetscInt i,idx=0; 142 const PetscScalar *xgen,*xnet; 143 PetscErrorCode ierr; 144 PetscScalar Efd,RF,VR,Vr,Vi,Vm; 145 146 PetscFunctionBegin; 147 148 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 149 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 150 151 ierr = VecGetArrayRead(Xgen,&xgen);CHKERRQ(ierr); 152 ierr = VecGetArrayRead(Xnet,&xnet);CHKERRQ(ierr); 153 154 /* Event for fault-on time */ 155 fvalue[0] = t - user->tfaulton; 156 /* Event for fault-off time */ 157 fvalue[1] = t - user->tfaultoff; 158 159 for (i=0; i < ngen; i++) { 160 Efd = xgen[idx+6]; 161 RF = xgen[idx+7]; 162 VR = xgen[idx+8]; 163 164 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 165 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 166 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 167 168 if (!VRatmax[i]) { 169 fvalue[2+2*i] = VRMAX[i] - VR; 170 } else { 171 fvalue[2+2*i] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; 172 } 173 if (!VRatmin[i]) { 174 fvalue[2+2*i+1] = VRMIN[i] - VR; 175 } else { 176 fvalue[2+2*i+1] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; 177 } 178 idx = idx+9; 179 } 180 ierr = VecRestoreArrayRead(Xgen,&xgen);CHKERRQ(ierr); 181 ierr = VecRestoreArrayRead(Xnet,&xnet);CHKERRQ(ierr); 182 183 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 184 185 PetscFunctionReturn(0); 186 } 187 188 PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec X,PetscBool forwardsolve,void* ctx) 189 { 190 Userctx *user=(Userctx*)ctx; 191 Vec Xgen,Xnet; 192 PetscScalar *xgen,*xnet; 193 PetscInt row_loc,col_loc; 194 PetscScalar val; 195 PetscErrorCode ierr; 196 PetscInt i,idx=0,event_num; 197 PetscScalar fvalue; 198 PetscScalar Efd, RF, VR; 199 PetscScalar Vr,Vi,Vm; 200 201 PetscFunctionBegin; 202 203 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 204 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 205 206 ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr); 207 ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr); 208 209 for (i=0; i < nevents; i++) { 210 if (event_list[i] == 0) { 211 /* Apply disturbance - resistive fault at user->faultbus */ 212 /* This is done by adding shunt conductance to the diagonal location 213 in the Ybus matrix */ 214 row_loc = 2*user->faultbus; col_loc = 2*user->faultbus+1; /* Location for G */ 215 val = 1/user->Rfault; 216 ierr = MatSetValues(user->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 217 row_loc = 2*user->faultbus+1; col_loc = 2*user->faultbus; /* Location for G */ 218 val = 1/user->Rfault; 219 ierr = MatSetValues(user->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 220 221 ierr = MatAssemblyBegin(user->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 222 ierr = MatAssemblyEnd(user->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 223 224 /* Solve the algebraic equations */ 225 ierr = SNESSolve(user->snes_alg,NULL,X);CHKERRQ(ierr); 226 } else if (event_list[i] == 1) { 227 /* Remove the fault */ 228 row_loc = 2*user->faultbus; col_loc = 2*user->faultbus+1; 229 val = -1/user->Rfault; 230 ierr = MatSetValues(user->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 231 row_loc = 2*user->faultbus+1; col_loc = 2*user->faultbus; 232 val = -1/user->Rfault; 233 ierr = MatSetValues(user->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);CHKERRQ(ierr); 234 235 ierr = MatAssemblyBegin(user->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 236 ierr = MatAssemblyEnd(user->Ybus,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 237 238 /* Solve the algebraic equations */ 239 ierr = SNESSolve(user->snes_alg,NULL,X);CHKERRQ(ierr); 240 241 /* Check the VR derivatives and reset flags if needed */ 242 for (i=0; i < ngen; i++) { 243 Efd = xgen[idx+6]; 244 RF = xgen[idx+7]; 245 VR = xgen[idx+8]; 246 247 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 248 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 249 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 250 251 if (VRatmax[i]) { 252 fvalue = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; 253 if (fvalue < 0) { 254 VRatmax[i] = 0; 255 ierr = PetscPrintf(PETSC_COMM_SELF,"VR[%d]: dVR_dt went negative on fault clearing at time %g\n",i,t);CHKERRQ(ierr); 256 } 257 } 258 if (VRatmin[i]) { 259 fvalue = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; 260 261 if (fvalue > 0) { 262 VRatmin[i] = 0; 263 ierr = PetscPrintf(PETSC_COMM_SELF,"VR[%d]: dVR_dt went positive on fault clearing at time %g\n",i,t);CHKERRQ(ierr); 264 } 265 } 266 idx = idx+9; 267 } 268 } else { 269 idx = (event_list[i]-2)/2; 270 event_num = (event_list[i]-2)%2; 271 if (event_num == 0) { /* Max VR */ 272 if (!VRatmax[idx]) { 273 VRatmax[idx] = 1; 274 ierr = PetscPrintf(PETSC_COMM_SELF,"VR[%d]: hit upper limit at time %g\n",idx,t);CHKERRQ(ierr); 275 } 276 else { 277 VRatmax[idx] = 0; 278 ierr = PetscPrintf(PETSC_COMM_SELF,"VR[%d]: freeing variable as dVR_dt is negative at time %g\n",idx,t);CHKERRQ(ierr); 279 } 280 } else { 281 if (!VRatmin[idx]) { 282 VRatmin[idx] = 1; 283 ierr = PetscPrintf(PETSC_COMM_SELF,"VR[%d]: hit lower limit at time %g\n",idx,t);CHKERRQ(ierr); 284 } 285 else { 286 VRatmin[idx] = 0; 287 ierr = PetscPrintf(PETSC_COMM_SELF,"VR[%d]: freeing variable as dVR_dt is positive at time %g\n",idx,t);CHKERRQ(ierr); 288 } 289 } 290 } 291 } 292 293 ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr); 294 ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr); 295 296 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 297 298 PetscFunctionReturn(0); 299 } 300 301 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */ 302 PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr, PetscScalar *Fi) 303 { 304 PetscFunctionBegin; 305 *Fr = Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta); 306 *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta); 307 PetscFunctionReturn(0); 308 } 309 310 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */ 311 PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd, PetscScalar *Fq) 312 { 313 PetscFunctionBegin; 314 *Fd = Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta); 315 *Fq = Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta); 316 PetscFunctionReturn(0); 317 } 318 319 /* Saves the solution at each time to a matrix */ 320 PetscErrorCode SaveSolution(TS ts) 321 { 322 PetscErrorCode ierr; 323 Userctx *user; 324 Vec X; 325 const PetscScalar *x; 326 PetscScalar *mat; 327 PetscInt idx; 328 PetscReal t; 329 330 PetscFunctionBegin; 331 ierr = TSGetApplicationContext(ts,&user);CHKERRQ(ierr); 332 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 333 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 334 idx = user->stepnum*(user->neqs_pgrid+1); 335 ierr = MatDenseGetArray(user->Sol,&mat);CHKERRQ(ierr); 336 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 337 mat[idx] = t; 338 ierr = PetscArraycpy(mat+idx+1,x,user->neqs_pgrid);CHKERRQ(ierr); 339 ierr = MatDenseRestoreArray(user->Sol,&mat);CHKERRQ(ierr); 340 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 341 user->stepnum++; 342 PetscFunctionReturn(0); 343 } 344 345 PetscErrorCode SetInitialGuess(Vec X,Userctx *user) 346 { 347 PetscErrorCode ierr; 348 Vec Xgen,Xnet; 349 PetscScalar *xgen; 350 const PetscScalar *xnet; 351 PetscInt i,idx=0; 352 PetscScalar Vr,Vi,IGr,IGi,Vm,Vm2; 353 PetscScalar Eqp,Edp,delta; 354 PetscScalar Efd,RF,VR; /* Exciter variables */ 355 PetscScalar Id,Iq; /* Generator dq axis currents */ 356 PetscScalar theta,Vd,Vq,SE; 357 358 PetscFunctionBegin; 359 M[0] = 2*H[0]/w_s; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s; 360 D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2]; 361 362 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 363 364 /* Network subsystem initialization */ 365 ierr = VecCopy(user->V0,Xnet);CHKERRQ(ierr); 366 367 /* Generator subsystem initialization */ 368 ierr = VecGetArrayWrite(Xgen,&xgen);CHKERRQ(ierr); 369 ierr = VecGetArrayRead(Xnet,&xnet);CHKERRQ(ierr); 370 371 for (i=0; i < ngen; i++) { 372 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 373 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 374 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 375 IGr = (Vr*PG[i] + Vi*QG[i])/Vm2; 376 IGi = (Vi*PG[i] - Vr*QG[i])/Vm2; 377 378 delta = PetscAtan2Real(Vi+Xq[i]*IGr,Vr-Xq[i]*IGi); /* Machine angle */ 379 380 theta = PETSC_PI/2.0 - delta; 381 382 Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); /* d-axis stator current */ 383 Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); /* q-axis stator current */ 384 385 Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta); 386 Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta); 387 388 Edp = Vd + Rs[i]*Id - Xqp[i]*Iq; /* d-axis transient EMF */ 389 Eqp = Vq + Rs[i]*Iq + Xdp[i]*Id; /* q-axis transient EMF */ 390 391 TM[i] = PG[i]; 392 393 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */ 394 xgen[idx] = Eqp; 395 xgen[idx+1] = Edp; 396 xgen[idx+2] = delta; 397 xgen[idx+3] = w_s; 398 399 idx = idx + 4; 400 401 xgen[idx] = Id; 402 xgen[idx+1] = Iq; 403 404 idx = idx + 2; 405 406 /* Exciter */ 407 Efd = Eqp + (Xd[i] - Xdp[i])*Id; 408 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 409 VR = KE[i]*Efd + SE; 410 RF = KF[i]*Efd/TF[i]; 411 412 xgen[idx] = Efd; 413 xgen[idx+1] = RF; 414 xgen[idx+2] = VR; 415 416 Vref[i] = Vm + (VR/KA[i]); 417 418 VRatmin[i] = VRatmax[i] = 0; 419 420 idx = idx + 3; 421 } 422 ierr = VecRestoreArrayWrite(Xgen,&xgen);CHKERRQ(ierr); 423 ierr = VecRestoreArrayRead(Xnet,&xnet);CHKERRQ(ierr); 424 425 /* ierr = VecView(Xgen,0);CHKERRQ(ierr); */ 426 ierr = DMCompositeGather(user->dmpgrid,INSERT_VALUES,X,Xgen,Xnet);CHKERRQ(ierr); 427 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 428 PetscFunctionReturn(0); 429 } 430 431 /* Computes F = [f(x,y);g(x,y)] */ 432 PetscErrorCode ResidualFunction(Vec X, Vec F, Userctx *user) 433 { 434 PetscErrorCode ierr; 435 Vec Xgen,Xnet,Fgen,Fnet; 436 const PetscScalar *xgen,*xnet; 437 PetscScalar *fgen,*fnet; 438 PetscInt i,idx=0; 439 PetscScalar Vr,Vi,Vm,Vm2; 440 PetscScalar Eqp,Edp,delta,w; /* Generator variables */ 441 PetscScalar Efd,RF,VR; /* Exciter variables */ 442 PetscScalar Id,Iq; /* Generator dq axis currents */ 443 PetscScalar Vd,Vq,SE; 444 PetscScalar IGr,IGi,IDr,IDi; 445 PetscScalar Zdq_inv[4],det; 446 PetscScalar PD,QD,Vm0,*v0; 447 PetscInt k; 448 449 PetscFunctionBegin; 450 ierr = VecZeroEntries(F);CHKERRQ(ierr); 451 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 452 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Fgen,&Fnet);CHKERRQ(ierr); 453 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 454 ierr = DMCompositeScatter(user->dmpgrid,F,Fgen,Fnet);CHKERRQ(ierr); 455 456 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here. 457 The generator current injection, IG, and load current injection, ID are added later 458 */ 459 /* Note that the values in Ybus are stored assuming the imaginary current balance 460 equation is ordered first followed by real current balance equation for each bus. 461 Thus imaginary current contribution goes in location 2*i, and 462 real current contribution in 2*i+1 463 */ 464 ierr = MatMult(user->Ybus,Xnet,Fnet);CHKERRQ(ierr); 465 466 ierr = VecGetArrayRead(Xgen,&xgen);CHKERRQ(ierr); 467 ierr = VecGetArrayRead(Xnet,&xnet);CHKERRQ(ierr); 468 ierr = VecGetArrayWrite(Fgen,&fgen);CHKERRQ(ierr); 469 ierr = VecGetArrayWrite(Fnet,&fnet);CHKERRQ(ierr); 470 471 /* Generator subsystem */ 472 for (i=0; i < ngen; i++) { 473 Eqp = xgen[idx]; 474 Edp = xgen[idx+1]; 475 delta = xgen[idx+2]; 476 w = xgen[idx+3]; 477 Id = xgen[idx+4]; 478 Iq = xgen[idx+5]; 479 Efd = xgen[idx+6]; 480 RF = xgen[idx+7]; 481 VR = xgen[idx+8]; 482 483 /* Generator differential equations */ 484 fgen[idx] = (-Eqp - (Xd[i] - Xdp[i])*Id + Efd)/Td0p[i]; 485 fgen[idx+1] = (-Edp + (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; 486 fgen[idx+2] = w - w_s; 487 fgen[idx+3] = (TM[i] - Edp*Id - Eqp*Iq - (Xqp[i] - Xdp[i])*Id*Iq - D[i]*(w - w_s))/M[i]; 488 489 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 490 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 491 492 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 493 /* Algebraic equations for stator currents */ 494 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 495 496 Zdq_inv[0] = Rs[i]/det; 497 Zdq_inv[1] = Xqp[i]/det; 498 Zdq_inv[2] = -Xdp[i]/det; 499 Zdq_inv[3] = Rs[i]/det; 500 501 fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 502 fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 503 504 /* Add generator current injection to network */ 505 ierr = dq2ri(Id,Iq,delta,&IGr,&IGi);CHKERRQ(ierr); 506 507 fnet[2*gbus[i]] -= IGi; 508 fnet[2*gbus[i]+1] -= IGr; 509 510 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 511 512 SE = k1[i]*PetscExpScalar(k2[i]*Efd); 513 514 /* Exciter differential equations */ 515 fgen[idx+6] = (-KE[i]*Efd - SE + VR)/TE[i]; 516 fgen[idx+7] = (-RF + KF[i]*Efd/TF[i])/TF[i]; 517 if (VRatmax[i]) fgen[idx+8] = VR - VRMAX[i]; 518 else if (VRatmin[i]) fgen[idx+8] = VRMIN[i] - VR; 519 else fgen[idx+8] = (-VR + KA[i]*RF - KA[i]*KF[i]*Efd/TF[i] + KA[i]*(Vref[i] - Vm))/TA[i]; 520 521 idx = idx + 9; 522 } 523 524 ierr = VecGetArray(user->V0,&v0);CHKERRQ(ierr); 525 for (i=0; i < nload; i++) { 526 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 527 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 528 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; 529 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 530 PD = QD = 0.0; 531 for (k=0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 532 for (k=0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 533 534 /* Load currents */ 535 IDr = (PD*Vr + QD*Vi)/Vm2; 536 IDi = (-QD*Vr + PD*Vi)/Vm2; 537 538 fnet[2*lbus[i]] += IDi; 539 fnet[2*lbus[i]+1] += IDr; 540 } 541 ierr = VecRestoreArray(user->V0,&v0);CHKERRQ(ierr); 542 543 ierr = VecRestoreArrayRead(Xgen,&xgen);CHKERRQ(ierr); 544 ierr = VecRestoreArrayRead(Xnet,&xnet);CHKERRQ(ierr); 545 ierr = VecRestoreArrayWrite(Fgen,&fgen);CHKERRQ(ierr); 546 ierr = VecRestoreArrayWrite(Fnet,&fnet);CHKERRQ(ierr); 547 548 ierr = DMCompositeGather(user->dmpgrid,INSERT_VALUES,F,Fgen,Fnet);CHKERRQ(ierr); 549 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 550 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Fgen,&Fnet);CHKERRQ(ierr); 551 PetscFunctionReturn(0); 552 } 553 554 /* f(x,y) 555 g(x,y) 556 */ 557 PetscErrorCode RHSFunction(TS ts,PetscReal t, Vec X, Vec F, void *ctx) 558 { 559 PetscErrorCode ierr; 560 Userctx *user=(Userctx*)ctx; 561 562 PetscFunctionBegin; 563 user->t = t; 564 ierr = ResidualFunction(X,F,user);CHKERRQ(ierr); 565 PetscFunctionReturn(0); 566 } 567 568 /* f(x,y) - \dot{x} 569 g(x,y) = 0 570 */ 571 PetscErrorCode IFunction(TS ts,PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx) 572 { 573 PetscErrorCode ierr; 574 PetscScalar *f; 575 const PetscScalar *xdot; 576 PetscInt i; 577 578 PetscFunctionBegin; 579 ierr = RHSFunction(ts,t,X,F,ctx);CHKERRQ(ierr); 580 ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 581 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 582 for (i=0;i < ngen;i++) { 583 f[9*i] -= xdot[9*i]; 584 f[9*i+1] -= xdot[9*i+1]; 585 f[9*i+2] -= xdot[9*i+2]; 586 f[9*i+3] -= xdot[9*i+3]; 587 f[9*i+6] -= xdot[9*i+6]; 588 f[9*i+7] -= xdot[9*i+7]; 589 f[9*i+8] -= xdot[9*i+8]; 590 } 591 ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 592 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 593 PetscFunctionReturn(0); 594 } 595 596 /* This function is used for solving the algebraic system only during fault on and 597 off times. It computes the entire F and then zeros out the part corresponding to 598 differential equations 599 F = [0;g(y)]; 600 */ 601 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx) 602 { 603 PetscErrorCode ierr; 604 Userctx *user=(Userctx*)ctx; 605 PetscScalar *f; 606 PetscInt i; 607 608 PetscFunctionBegin; 609 ierr = ResidualFunction(X,F,user);CHKERRQ(ierr); 610 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 611 for (i=0; i < ngen; i++) { 612 f[9*i] = 0; 613 f[9*i+1] = 0; 614 f[9*i+2] = 0; 615 f[9*i+3] = 0; 616 f[9*i+6] = 0; 617 f[9*i+7] = 0; 618 f[9*i+8] = 0; 619 } 620 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 621 PetscFunctionReturn(0); 622 } 623 624 PetscErrorCode PostStage(TS ts, PetscReal t, PetscInt i, Vec *X) 625 { 626 PetscErrorCode ierr; 627 Userctx *user; 628 629 PetscFunctionBegin; 630 ierr = TSGetApplicationContext(ts,&user);CHKERRQ(ierr); 631 ierr = SNESSolve(user->snes_alg,NULL,X[i]);CHKERRQ(ierr); 632 PetscFunctionReturn(0); 633 } 634 635 PetscErrorCode PostEvaluate(TS ts) 636 { 637 PetscErrorCode ierr; 638 Userctx *user; 639 Vec X; 640 641 PetscFunctionBegin; 642 ierr = TSGetApplicationContext(ts,&user);CHKERRQ(ierr); 643 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 644 ierr = SNESSolve(user->snes_alg,NULL,X);CHKERRQ(ierr); 645 PetscFunctionReturn(0); 646 } 647 648 649 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user) 650 { 651 PetscErrorCode ierr; 652 PetscInt *d_nnz; 653 PetscInt i,idx=0,start=0; 654 PetscInt ncols; 655 656 PetscFunctionBegin; 657 ierr = PetscMalloc1(user->neqs_pgrid,&d_nnz);CHKERRQ(ierr); 658 for (i=0; i<user->neqs_pgrid; i++) d_nnz[i] = 0; 659 /* Generator subsystem */ 660 for (i=0; i < ngen; i++) { 661 662 d_nnz[idx] += 3; 663 d_nnz[idx+1] += 2; 664 d_nnz[idx+2] += 2; 665 d_nnz[idx+3] += 5; 666 d_nnz[idx+4] += 6; 667 d_nnz[idx+5] += 6; 668 669 d_nnz[user->neqs_gen+2*gbus[i]] += 3; 670 d_nnz[user->neqs_gen+2*gbus[i]+1] += 3; 671 672 d_nnz[idx+6] += 2; 673 d_nnz[idx+7] += 2; 674 d_nnz[idx+8] += 5; 675 676 idx = idx + 9; 677 } 678 start = user->neqs_gen; 679 680 for (i=0; i < nbus; i++) { 681 ierr = MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 682 d_nnz[start+2*i] += ncols; 683 d_nnz[start+2*i+1] += ncols; 684 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 685 } 686 ierr = MatSeqAIJSetPreallocation(J,0,d_nnz);CHKERRQ(ierr); 687 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 688 PetscFunctionReturn(0); 689 } 690 691 /* 692 J = [df_dx, df_dy 693 dg_dx, dg_dy] 694 */ 695 PetscErrorCode ResidualJacobian(Vec X,Mat J,Mat B,void *ctx) 696 { 697 PetscErrorCode ierr; 698 Userctx *user = (Userctx*)ctx; 699 Vec Xgen,Xnet; 700 const PetscScalar *xgen,*xnet; 701 PetscInt i,idx=0; 702 PetscScalar Vr,Vi,Vm,Vm2; 703 PetscScalar Eqp,Edp,delta; /* Generator variables */ 704 PetscScalar Efd; 705 PetscScalar Id,Iq; /* Generator dq axis currents */ 706 PetscScalar Vd,Vq; 707 PetscScalar val[10]; 708 PetscInt row[2],col[10]; 709 PetscInt net_start = user->neqs_gen; 710 PetscScalar Zdq_inv[4],det; 711 PetscScalar dVd_dVr,dVd_dVi,dVq_dVr,dVq_dVi,dVd_ddelta,dVq_ddelta; 712 PetscScalar dIGr_ddelta,dIGi_ddelta,dIGr_dId,dIGr_dIq,dIGi_dId,dIGi_dIq; 713 PetscScalar dSE_dEfd; 714 PetscScalar dVm_dVd,dVm_dVq,dVm_dVr,dVm_dVi; 715 PetscInt ncols; 716 const PetscInt *cols; 717 const PetscScalar *yvals; 718 PetscInt k; 719 PetscScalar PD,QD,Vm0,Vm4; 720 const PetscScalar *v0; 721 PetscScalar dPD_dVr,dPD_dVi,dQD_dVr,dQD_dVi; 722 PetscScalar dIDr_dVr,dIDr_dVi,dIDi_dVr,dIDi_dVi; 723 724 725 PetscFunctionBegin; 726 ierr = MatZeroEntries(B);CHKERRQ(ierr); 727 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 728 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 729 730 ierr = VecGetArrayRead(Xgen,&xgen);CHKERRQ(ierr); 731 ierr = VecGetArrayRead(Xnet,&xnet);CHKERRQ(ierr); 732 733 /* Generator subsystem */ 734 for (i=0; i < ngen; i++) { 735 Eqp = xgen[idx]; 736 Edp = xgen[idx+1]; 737 delta = xgen[idx+2]; 738 Id = xgen[idx+4]; 739 Iq = xgen[idx+5]; 740 Efd = xgen[idx+6]; 741 742 /* fgen[idx] = (-Eqp - (Xd[i] - Xdp[i])*Id + Efd)/Td0p[i]; */ 743 row[0] = idx; 744 col[0] = idx; col[1] = idx+4; col[2] = idx+6; 745 val[0] = -1/ Td0p[i]; val[1] = -(Xd[i] - Xdp[i])/ Td0p[i]; val[2] = 1/Td0p[i]; 746 747 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 748 749 /* fgen[idx+1] = (-Edp + (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */ 750 row[0] = idx + 1; 751 col[0] = idx + 1; col[1] = idx+5; 752 val[0] = -1/Tq0p[i]; val[1] = (Xq[i] - Xqp[i])/Tq0p[i]; 753 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 754 755 /* fgen[idx+2] = w - w_s; */ 756 row[0] = idx + 2; 757 col[0] = idx + 2; col[1] = idx + 3; 758 val[0] = 0; val[1] = 1; 759 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 760 761 /* fgen[idx+3] = (TM[i] - Edp*Id - Eqp*Iq - (Xqp[i] - Xdp[i])*Id*Iq - D[i]*(w - w_s))/M[i]; */ 762 row[0] = idx + 3; 763 col[0] = idx; col[1] = idx + 1; col[2] = idx + 3; col[3] = idx + 4; col[4] = idx + 5; 764 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]; 765 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 766 767 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 768 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 769 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 770 771 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 772 773 Zdq_inv[0] = Rs[i]/det; 774 Zdq_inv[1] = Xqp[i]/det; 775 Zdq_inv[2] = -Xdp[i]/det; 776 Zdq_inv[3] = Rs[i]/det; 777 778 dVd_dVr = PetscSinScalar(delta); dVd_dVi = -PetscCosScalar(delta); 779 dVq_dVr = PetscCosScalar(delta); dVq_dVi = PetscSinScalar(delta); 780 dVd_ddelta = Vr*PetscCosScalar(delta) + Vi*PetscSinScalar(delta); 781 dVq_ddelta = -Vr*PetscSinScalar(delta) + Vi*PetscCosScalar(delta); 782 783 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */ 784 row[0] = idx+4; 785 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 786 val[0] = -Zdq_inv[1]; val[1] = -Zdq_inv[0]; val[2] = Zdq_inv[0]*dVd_ddelta + Zdq_inv[1]*dVq_ddelta; 787 col[3] = idx + 4; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 788 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; 789 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 790 791 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */ 792 row[0] = idx+5; 793 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 794 val[0] = -Zdq_inv[3]; val[1] = -Zdq_inv[2]; val[2] = Zdq_inv[2]*dVd_ddelta + Zdq_inv[3]*dVq_ddelta; 795 col[3] = idx + 5; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 796 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; 797 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 798 799 dIGr_ddelta = Id*PetscCosScalar(delta) - Iq*PetscSinScalar(delta); 800 dIGi_ddelta = Id*PetscSinScalar(delta) + Iq*PetscCosScalar(delta); 801 dIGr_dId = PetscSinScalar(delta); dIGr_dIq = PetscCosScalar(delta); 802 dIGi_dId = -PetscCosScalar(delta); dIGi_dIq = PetscSinScalar(delta); 803 804 /* fnet[2*gbus[i]] -= IGi; */ 805 row[0] = net_start + 2*gbus[i]; 806 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 807 val[0] = -dIGi_ddelta; val[1] = -dIGi_dId; val[2] = -dIGi_dIq; 808 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 809 810 /* fnet[2*gbus[i]+1] -= IGr; */ 811 row[0] = net_start + 2*gbus[i]+1; 812 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 813 val[0] = -dIGr_ddelta; val[1] = -dIGr_dId; val[2] = -dIGr_dIq; 814 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 815 816 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 817 818 /* fgen[idx+6] = (-KE[i]*Efd - SE + VR)/TE[i]; */ 819 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */ 820 821 dSE_dEfd = k1[i]*k2[i]*PetscExpScalar(k2[i]*Efd); 822 823 row[0] = idx + 6; 824 col[0] = idx + 6; col[1] = idx + 8; 825 val[0] = (-KE[i] - dSE_dEfd)/TE[i]; val[1] = 1/TE[i]; 826 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 827 828 /* Exciter differential equations */ 829 830 /* fgen[idx+7] = (-RF + KF[i]*Efd/TF[i])/TF[i]; */ 831 row[0] = idx + 7; 832 col[0] = idx + 6; col[1] = idx + 7; 833 val[0] = (KF[i]/TF[i])/TF[i]; val[1] = -1/TF[i]; 834 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 835 836 /* fgen[idx+8] = (-VR + KA[i]*RF - KA[i]*KF[i]*Efd/TF[i] + KA[i]*(Vref[i] - Vm))/TA[i]; */ 837 /* Vm = (Vd^2 + Vq^2)^0.5; */ 838 839 row[0] = idx + 8; 840 if (VRatmax[i]) { 841 col[0] = idx + 8; val[0] = 1.0; 842 ierr = MatSetValues(J,1,row,1,col,val,INSERT_VALUES);CHKERRQ(ierr); 843 } else if (VRatmin[i]) { 844 col[0] = idx + 8; val[0] = -1.0; 845 ierr = MatSetValues(J,1,row,1,col,val,INSERT_VALUES);CHKERRQ(ierr); 846 } else { 847 dVm_dVd = Vd/Vm; dVm_dVq = Vq/Vm; 848 dVm_dVr = dVm_dVd*dVd_dVr + dVm_dVq*dVq_dVr; 849 dVm_dVi = dVm_dVd*dVd_dVi + dVm_dVq*dVq_dVi; 850 row[0] = idx + 8; 851 col[0] = idx + 6; col[1] = idx + 7; col[2] = idx + 8; 852 val[0] = -(KA[i]*KF[i]/TF[i])/TA[i]; val[1] = KA[i]/TA[i]; val[2] = -1/TA[i]; 853 col[3] = net_start + 2*gbus[i]; col[4] = net_start + 2*gbus[i]+1; 854 val[3] = -KA[i]*dVm_dVr/TA[i]; val[4] = -KA[i]*dVm_dVi/TA[i]; 855 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 856 } 857 idx = idx + 9; 858 } 859 860 for (i=0; i<nbus; i++) { 861 ierr = MatGetRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 862 row[0] = net_start + 2*i; 863 for (k=0; k<ncols; k++) { 864 col[k] = net_start + cols[k]; 865 val[k] = yvals[k]; 866 } 867 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 868 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 869 870 ierr = MatGetRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 871 row[0] = net_start + 2*i+1; 872 for (k=0; k<ncols; k++) { 873 col[k] = net_start + cols[k]; 874 val[k] = yvals[k]; 875 } 876 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 877 ierr = MatRestoreRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 878 } 879 880 ierr = MatAssemblyBegin(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 881 ierr = MatAssemblyEnd(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 882 883 ierr = VecGetArrayRead(user->V0,&v0);CHKERRQ(ierr); 884 for (i=0; i < nload; i++) { 885 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 886 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 887 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; Vm4 = Vm2*Vm2; 888 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 889 PD = QD = 0.0; 890 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0; 891 for (k=0; k < ld_nsegsp[i]; k++) { 892 PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 893 dPD_dVr += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vr*PetscPowScalar(Vm,(ld_betap[k]-2)); 894 dPD_dVi += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vi*PetscPowScalar(Vm,(ld_betap[k]-2)); 895 } 896 for (k=0; k < ld_nsegsq[i]; k++) { 897 QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 898 dQD_dVr += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vr*PetscPowScalar(Vm,(ld_betaq[k]-2)); 899 dQD_dVi += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vi*PetscPowScalar(Vm,(ld_betaq[k]-2)); 900 } 901 902 /* IDr = (PD*Vr + QD*Vi)/Vm2; */ 903 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */ 904 905 dIDr_dVr = (dPD_dVr*Vr + dQD_dVr*Vi + PD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vr)/Vm4; 906 dIDr_dVi = (dPD_dVi*Vr + dQD_dVi*Vi + QD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vi)/Vm4; 907 908 dIDi_dVr = (-dQD_dVr*Vr + dPD_dVr*Vi - QD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vr)/Vm4; 909 dIDi_dVi = (-dQD_dVi*Vr + dPD_dVi*Vi + PD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vi)/Vm4; 910 911 912 /* fnet[2*lbus[i]] += IDi; */ 913 row[0] = net_start + 2*lbus[i]; 914 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 915 val[0] = dIDi_dVr; val[1] = dIDi_dVi; 916 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 917 /* fnet[2*lbus[i]+1] += IDr; */ 918 row[0] = net_start + 2*lbus[i]+1; 919 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 920 val[0] = dIDr_dVr; val[1] = dIDr_dVi; 921 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 922 } 923 ierr = VecRestoreArrayRead(user->V0,&v0);CHKERRQ(ierr); 924 925 ierr = VecRestoreArrayRead(Xgen,&xgen);CHKERRQ(ierr); 926 ierr = VecRestoreArrayRead(Xnet,&xnet);CHKERRQ(ierr); 927 928 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 929 930 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 931 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 932 PetscFunctionReturn(0); 933 } 934 935 /* 936 J = [I, 0 937 dg_dx, dg_dy] 938 */ 939 PetscErrorCode AlgJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx) 940 { 941 PetscErrorCode ierr; 942 Userctx *user=(Userctx*)ctx; 943 944 PetscFunctionBegin; 945 ierr = ResidualJacobian(X,A,B,ctx);CHKERRQ(ierr); 946 ierr = MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);CHKERRQ(ierr); 947 ierr = MatZeroRowsIS(A,user->is_diff,1.0,NULL,NULL);CHKERRQ(ierr); 948 PetscFunctionReturn(0); 949 } 950 951 /* 952 J = [-df_dx, -df_dy 953 dg_dx, dg_dy] 954 */ 955 956 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec X,Mat A,Mat B,void *ctx) 957 { 958 PetscErrorCode ierr; 959 Userctx *user=(Userctx*)ctx; 960 961 PetscFunctionBegin; 962 user->t = t; 963 964 ierr = ResidualJacobian(X,A,B,user);CHKERRQ(ierr); 965 966 PetscFunctionReturn(0); 967 } 968 969 /* 970 J = [df_dx-aI, df_dy 971 dg_dx, dg_dy] 972 */ 973 974 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx *user) 975 { 976 PetscErrorCode ierr; 977 PetscScalar atmp = (PetscScalar) a; 978 PetscInt i,row; 979 980 PetscFunctionBegin; 981 user->t = t; 982 atmp *= -1; 983 984 ierr = RHSJacobian(ts,t,X,A,B,user);CHKERRQ(ierr); 985 for (i=0;i < ngen;i++) { 986 row = 9*i; 987 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 988 row = 9*i+1; 989 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 990 row = 9*i+2; 991 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 992 row = 9*i+3; 993 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 994 row = 9*i+6; 995 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 996 row = 9*i+7; 997 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 998 row = 9*i+8; 999 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 1000 } 1001 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1002 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1003 PetscFunctionReturn(0); 1004 } 1005 1006 int main(int argc,char **argv) 1007 { 1008 TS ts; 1009 SNES snes_alg; 1010 PetscErrorCode ierr; 1011 PetscMPIInt size; 1012 Userctx user; 1013 PetscViewer Xview,Ybusview,viewer; 1014 Vec X,F_alg; 1015 Mat J,A; 1016 PetscInt i,idx,*idx2; 1017 Vec Xdot; 1018 PetscScalar *x,*mat,*amat; 1019 const PetscScalar *rmat; 1020 Vec vatol; 1021 PetscInt *direction; 1022 PetscBool *terminate; 1023 const PetscInt *idx3; 1024 PetscScalar *vatoli; 1025 PetscInt k; 1026 1027 1028 ierr = PetscInitialize(&argc,&argv,"petscoptions",help);if (ierr) return ierr; 1029 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 1030 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 1031 1032 user.neqs_gen = 9*ngen; /* # eqs. for generator subsystem */ 1033 user.neqs_net = 2*nbus; /* # eqs. for network subsystem */ 1034 user.neqs_pgrid = user.neqs_gen + user.neqs_net; 1035 1036 /* Create indices for differential and algebraic equations */ 1037 1038 ierr = PetscMalloc1(7*ngen,&idx2);CHKERRQ(ierr); 1039 for (i=0; i<ngen; i++) { 1040 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; 1041 idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8; 1042 } 1043 ierr = ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff);CHKERRQ(ierr); 1044 ierr = ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg);CHKERRQ(ierr); 1045 ierr = PetscFree(idx2);CHKERRQ(ierr); 1046 1047 /* Read initial voltage vector and Ybus */ 1048 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"X.bin",FILE_MODE_READ,&Xview);CHKERRQ(ierr); 1049 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ybus.bin",FILE_MODE_READ,&Ybusview);CHKERRQ(ierr); 1050 1051 ierr = VecCreate(PETSC_COMM_WORLD,&user.V0);CHKERRQ(ierr); 1052 ierr = VecSetSizes(user.V0,PETSC_DECIDE,user.neqs_net);CHKERRQ(ierr); 1053 ierr = VecLoad(user.V0,Xview);CHKERRQ(ierr); 1054 1055 ierr = MatCreate(PETSC_COMM_WORLD,&user.Ybus);CHKERRQ(ierr); 1056 ierr = MatSetSizes(user.Ybus,PETSC_DECIDE,PETSC_DECIDE,user.neqs_net,user.neqs_net);CHKERRQ(ierr); 1057 ierr = MatSetType(user.Ybus,MATBAIJ);CHKERRQ(ierr); 1058 /* ierr = MatSetBlockSize(user.Ybus,2);CHKERRQ(ierr); */ 1059 ierr = MatLoad(user.Ybus,Ybusview);CHKERRQ(ierr); 1060 1061 /* Set run time options */ 1062 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");CHKERRQ(ierr); 1063 { 1064 user.tfaulton = 1.0; 1065 user.tfaultoff = 1.2; 1066 user.Rfault = 0.0001; 1067 user.setisdiff = PETSC_FALSE; 1068 user.semiexplicit = PETSC_FALSE; 1069 user.faultbus = 8; 1070 ierr = PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL);CHKERRQ(ierr); 1071 ierr = PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL);CHKERRQ(ierr); 1072 ierr = PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL);CHKERRQ(ierr); 1073 user.t0 = 0.0; 1074 user.tmax = 5.0; 1075 ierr = PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL);CHKERRQ(ierr); 1076 ierr = PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL);CHKERRQ(ierr); 1077 ierr = PetscOptionsBool("-setisdiff","","",user.setisdiff,&user.setisdiff,NULL);CHKERRQ(ierr); 1078 ierr = PetscOptionsBool("-dae_semiexplicit","","",user.semiexplicit,&user.semiexplicit,NULL);CHKERRQ(ierr); 1079 } 1080 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1081 1082 ierr = PetscViewerDestroy(&Xview);CHKERRQ(ierr); 1083 ierr = PetscViewerDestroy(&Ybusview);CHKERRQ(ierr); 1084 1085 /* Create DMs for generator and network subsystems */ 1086 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); 1087 ierr = DMSetOptionsPrefix(user.dmgen,"dmgen_");CHKERRQ(ierr); 1088 ierr = DMSetFromOptions(user.dmgen);CHKERRQ(ierr); 1089 ierr = DMSetUp(user.dmgen);CHKERRQ(ierr); 1090 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); 1091 ierr = DMSetOptionsPrefix(user.dmnet,"dmnet_");CHKERRQ(ierr); 1092 ierr = DMSetFromOptions(user.dmnet);CHKERRQ(ierr); 1093 ierr = DMSetUp(user.dmnet);CHKERRQ(ierr); 1094 /* Create a composite DM packer and add the two DMs */ 1095 ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid);CHKERRQ(ierr); 1096 ierr = DMSetOptionsPrefix(user.dmpgrid,"pgrid_");CHKERRQ(ierr); 1097 ierr = DMCompositeAddDM(user.dmpgrid,user.dmgen);CHKERRQ(ierr); 1098 ierr = DMCompositeAddDM(user.dmpgrid,user.dmnet);CHKERRQ(ierr); 1099 1100 ierr = DMCreateGlobalVector(user.dmpgrid,&X);CHKERRQ(ierr); 1101 1102 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 1103 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,user.neqs_pgrid,user.neqs_pgrid);CHKERRQ(ierr); 1104 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 1105 ierr = PreallocateJacobian(J,&user);CHKERRQ(ierr); 1106 1107 /* Create matrix to save solutions at each time step */ 1108 user.stepnum = 0; 1109 1110 ierr = MatCreateSeqDense(PETSC_COMM_SELF,user.neqs_pgrid+1,1002,NULL,&user.Sol);CHKERRQ(ierr); 1111 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1112 Create timestepping solver context 1113 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1114 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 1115 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 1116 if (user.semiexplicit) { 1117 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr); 1118 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); 1119 ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,&user);CHKERRQ(ierr); 1120 } else { 1121 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 1122 ierr = TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);CHKERRQ(ierr); 1123 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&user);CHKERRQ(ierr); 1124 ierr = TSSetIJacobian(ts,J,J,(TSIJacobian)IJacobian,&user);CHKERRQ(ierr); 1125 } 1126 ierr = TSSetApplicationContext(ts,&user);CHKERRQ(ierr); 1127 1128 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1129 Set initial conditions 1130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1131 ierr = SetInitialGuess(X,&user);CHKERRQ(ierr); 1132 /* Just to set up the Jacobian structure */ 1133 1134 ierr = VecDuplicate(X,&Xdot);CHKERRQ(ierr); 1135 ierr = IJacobian(ts,0.0,X,Xdot,0.0,J,J,&user);CHKERRQ(ierr); 1136 ierr = VecDestroy(&Xdot);CHKERRQ(ierr); 1137 1138 /* Save initial solution */ 1139 1140 idx=user.stepnum*(user.neqs_pgrid+1); 1141 ierr = MatDenseGetArray(user.Sol,&mat);CHKERRQ(ierr); 1142 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 1143 1144 mat[idx] = 0.0; 1145 1146 ierr = PetscArraycpy(mat+idx+1,x,user.neqs_pgrid);CHKERRQ(ierr); 1147 ierr = MatDenseRestoreArray(user.Sol,&mat);CHKERRQ(ierr); 1148 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 1149 user.stepnum++; 1150 1151 ierr = TSSetMaxTime(ts,user.tmax);CHKERRQ(ierr); 1152 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 1153 ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr); 1154 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 1155 ierr = TSSetPostStep(ts,SaveSolution);CHKERRQ(ierr); 1156 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 1157 1158 ierr = PetscMalloc1((2*ngen+2),&direction);CHKERRQ(ierr); 1159 ierr = PetscMalloc1((2*ngen+2),&terminate);CHKERRQ(ierr); 1160 direction[0] = direction[1] = 1; 1161 terminate[0] = terminate[1] = PETSC_FALSE; 1162 for (i=0; i < ngen;i++) { 1163 direction[2+2*i] = -1; direction[2+2*i+1] = 1; 1164 terminate[2+2*i] = terminate[2+2*i+1] = PETSC_FALSE; 1165 } 1166 1167 ierr = TSSetEventHandler(ts,2*ngen+2,direction,terminate,EventFunction,PostEventFunction,(void*)&user);CHKERRQ(ierr); 1168 1169 if (user.semiexplicit) { 1170 /* Use a semi-explicit approach with the time-stepping done by an explicit method and the 1171 algrebraic part solved via PostStage and PostEvaluate callbacks 1172 */ 1173 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr); 1174 ierr = TSSetPostStage(ts,PostStage);CHKERRQ(ierr); 1175 ierr = TSSetPostEvaluate(ts,PostEvaluate);CHKERRQ(ierr); 1176 } 1177 1178 1179 if (user.setisdiff) { 1180 /* Create vector of absolute tolerances and set the algebraic part to infinity */ 1181 ierr = VecDuplicate(X,&vatol);CHKERRQ(ierr); 1182 ierr = VecSet(vatol,100000.0);CHKERRQ(ierr); 1183 ierr = VecGetArray(vatol,&vatoli);CHKERRQ(ierr); 1184 ierr = ISGetIndices(user.is_diff,&idx3);CHKERRQ(ierr); 1185 for (k=0; k < 7*ngen; k++) vatoli[idx3[k]] = 1e-2; 1186 ierr = VecRestoreArray(vatol,&vatoli);CHKERRQ(ierr); 1187 } 1188 1189 /* Create the nonlinear solver for solving the algebraic system */ 1190 /* Note that although the algebraic system needs to be solved only for 1191 Idq and V, we reuse the entire system including xgen. The xgen 1192 variables are held constant by setting their residuals to 0 and 1193 putting a 1 on the Jacobian diagonal for xgen rows 1194 */ 1195 1196 ierr = VecDuplicate(X,&F_alg);CHKERRQ(ierr); 1197 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_alg);CHKERRQ(ierr); 1198 ierr = SNESSetFunction(snes_alg,F_alg,AlgFunction,&user);CHKERRQ(ierr); 1199 ierr = SNESSetJacobian(snes_alg,J,J,AlgJacobian,&user);CHKERRQ(ierr); 1200 ierr = SNESSetFromOptions(snes_alg);CHKERRQ(ierr); 1201 1202 user.snes_alg=snes_alg; 1203 1204 /* Solve */ 1205 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1206 1207 ierr = MatAssemblyBegin(user.Sol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1208 ierr = MatAssemblyEnd(user.Sol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1209 1210 ierr = MatCreateSeqDense(PETSC_COMM_SELF,user.neqs_pgrid+1,user.stepnum,NULL,&A);CHKERRQ(ierr); 1211 ierr = MatDenseGetArrayRead(user.Sol,&rmat);CHKERRQ(ierr); 1212 ierr = MatDenseGetArray(A,&amat);CHKERRQ(ierr); 1213 ierr = PetscArraycpy(amat,rmat,user.stepnum*(user.neqs_pgrid+1));CHKERRQ(ierr); 1214 ierr = MatDenseRestoreArray(A,&amat);CHKERRQ(ierr); 1215 ierr = MatDenseRestoreArrayRead(user.Sol,&rmat);CHKERRQ(ierr); 1216 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"out.bin",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 1217 ierr = MatView(A,viewer);CHKERRQ(ierr); 1218 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1219 ierr = MatDestroy(&A);CHKERRQ(ierr); 1220 1221 ierr = PetscFree(direction);CHKERRQ(ierr); 1222 ierr = PetscFree(terminate);CHKERRQ(ierr); 1223 ierr = SNESDestroy(&snes_alg);CHKERRQ(ierr); 1224 ierr = VecDestroy(&F_alg);CHKERRQ(ierr); 1225 ierr = MatDestroy(&J);CHKERRQ(ierr); 1226 ierr = MatDestroy(&user.Ybus);CHKERRQ(ierr); 1227 ierr = MatDestroy(&user.Sol);CHKERRQ(ierr); 1228 ierr = VecDestroy(&X);CHKERRQ(ierr); 1229 ierr = VecDestroy(&user.V0);CHKERRQ(ierr); 1230 ierr = DMDestroy(&user.dmgen);CHKERRQ(ierr); 1231 ierr = DMDestroy(&user.dmnet);CHKERRQ(ierr); 1232 ierr = DMDestroy(&user.dmpgrid);CHKERRQ(ierr); 1233 ierr = ISDestroy(&user.is_diff);CHKERRQ(ierr); 1234 ierr = ISDestroy(&user.is_alg);CHKERRQ(ierr); 1235 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1236 if (user.setisdiff) { 1237 ierr = VecDestroy(&vatol);CHKERRQ(ierr); 1238 } 1239 ierr = PetscFinalize(); 1240 return ierr; 1241 } 1242 1243 /*TEST 1244 1245 build: 1246 requires: double !complex !define(PETSC_USE_64BIT_INDICES) 1247 1248 test: 1249 suffix: implicit 1250 args: -ts_monitor -snes_monitor_short 1251 localrunfiles: petscoptions X.bin Ybus.bin 1252 1253 test: 1254 suffix: semiexplicit 1255 args: -ts_monitor -snes_monitor_short -dae_semiexplicit -ts_rk_type 2a 1256 localrunfiles: petscoptions X.bin Ybus.bin 1257 1258 test: 1259 suffix: steprestart 1260 args: -ts_monitor -snes_monitor_short -ts_type arkimex 1261 localrunfiles: petscoptions X.bin Ybus.bin 1262 1263 TEST*/ 1264