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 = VecScale(F,-1.0);CHKERRQ(ierr); 581 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 582 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 583 for (i=0;i < ngen;i++) { 584 f[9*i] += xdot[9*i]; 585 f[9*i+1] += xdot[9*i+1]; 586 f[9*i+2] += xdot[9*i+2]; 587 f[9*i+3] += xdot[9*i+3]; 588 f[9*i+6] += xdot[9*i+6]; 589 f[9*i+7] += xdot[9*i+7]; 590 f[9*i+8] += xdot[9*i+8]; 591 } 592 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 593 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 594 PetscFunctionReturn(0); 595 } 596 597 /* This function is used for solving the algebraic system only during fault on and 598 off times. It computes the entire F and then zeros out the part corresponding to 599 differential equations 600 F = [0;g(y)]; 601 */ 602 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx) 603 { 604 PetscErrorCode ierr; 605 Userctx *user=(Userctx*)ctx; 606 PetscScalar *f; 607 PetscInt i; 608 609 PetscFunctionBegin; 610 ierr = ResidualFunction(X,F,user);CHKERRQ(ierr); 611 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 612 for (i=0; i < ngen; i++) { 613 f[9*i] = 0; 614 f[9*i+1] = 0; 615 f[9*i+2] = 0; 616 f[9*i+3] = 0; 617 f[9*i+6] = 0; 618 f[9*i+7] = 0; 619 f[9*i+8] = 0; 620 } 621 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode PostStage(TS ts, PetscReal t, PetscInt i, Vec *X) 626 { 627 PetscErrorCode ierr; 628 Userctx *user; 629 630 PetscFunctionBegin; 631 ierr = TSGetApplicationContext(ts,&user);CHKERRQ(ierr); 632 ierr = SNESSolve(user->snes_alg,NULL,X[i]);CHKERRQ(ierr); 633 PetscFunctionReturn(0); 634 } 635 636 PetscErrorCode PostEvaluate(TS ts) 637 { 638 PetscErrorCode ierr; 639 Userctx *user; 640 Vec X; 641 642 PetscFunctionBegin; 643 ierr = TSGetApplicationContext(ts,&user);CHKERRQ(ierr); 644 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 645 ierr = SNESSolve(user->snes_alg,NULL,X);CHKERRQ(ierr); 646 PetscFunctionReturn(0); 647 } 648 649 650 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user) 651 { 652 PetscErrorCode ierr; 653 PetscInt *d_nnz; 654 PetscInt i,idx=0,start=0; 655 PetscInt ncols; 656 657 PetscFunctionBegin; 658 ierr = PetscMalloc1(user->neqs_pgrid,&d_nnz);CHKERRQ(ierr); 659 for (i=0; i<user->neqs_pgrid; i++) d_nnz[i] = 0; 660 /* Generator subsystem */ 661 for (i=0; i < ngen; i++) { 662 663 d_nnz[idx] += 3; 664 d_nnz[idx+1] += 2; 665 d_nnz[idx+2] += 2; 666 d_nnz[idx+3] += 5; 667 d_nnz[idx+4] += 6; 668 d_nnz[idx+5] += 6; 669 670 d_nnz[user->neqs_gen+2*gbus[i]] += 3; 671 d_nnz[user->neqs_gen+2*gbus[i]+1] += 3; 672 673 d_nnz[idx+6] += 2; 674 d_nnz[idx+7] += 2; 675 d_nnz[idx+8] += 5; 676 677 idx = idx + 9; 678 } 679 start = user->neqs_gen; 680 681 for (i=0; i < nbus; i++) { 682 ierr = MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 683 d_nnz[start+2*i] += ncols; 684 d_nnz[start+2*i+1] += ncols; 685 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr); 686 } 687 ierr = MatSeqAIJSetPreallocation(J,0,d_nnz);CHKERRQ(ierr); 688 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 689 PetscFunctionReturn(0); 690 } 691 692 /* 693 J = [df_dx, df_dy 694 dg_dx, dg_dy] 695 */ 696 PetscErrorCode ResidualJacobian(Vec X,Mat J,Mat B,void *ctx) 697 { 698 PetscErrorCode ierr; 699 Userctx *user = (Userctx*)ctx; 700 Vec Xgen,Xnet; 701 const PetscScalar *xgen,*xnet; 702 PetscInt i,idx=0; 703 PetscScalar Vr,Vi,Vm,Vm2; 704 PetscScalar Eqp,Edp,delta; /* Generator variables */ 705 PetscScalar Efd; 706 PetscScalar Id,Iq; /* Generator dq axis currents */ 707 PetscScalar Vd,Vq; 708 PetscScalar val[10]; 709 PetscInt row[2],col[10]; 710 PetscInt net_start = user->neqs_gen; 711 PetscScalar Zdq_inv[4],det; 712 PetscScalar dVd_dVr,dVd_dVi,dVq_dVr,dVq_dVi,dVd_ddelta,dVq_ddelta; 713 PetscScalar dIGr_ddelta,dIGi_ddelta,dIGr_dId,dIGr_dIq,dIGi_dId,dIGi_dIq; 714 PetscScalar dSE_dEfd; 715 PetscScalar dVm_dVd,dVm_dVq,dVm_dVr,dVm_dVi; 716 PetscInt ncols; 717 const PetscInt *cols; 718 const PetscScalar *yvals; 719 PetscInt k; 720 PetscScalar PD,QD,Vm0,Vm4; 721 const PetscScalar *v0; 722 PetscScalar dPD_dVr,dPD_dVi,dQD_dVr,dQD_dVi; 723 PetscScalar dIDr_dVr,dIDr_dVi,dIDi_dVr,dIDi_dVi; 724 725 726 PetscFunctionBegin; 727 ierr = MatZeroEntries(B);CHKERRQ(ierr); 728 ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 729 ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr); 730 731 ierr = VecGetArrayRead(Xgen,&xgen);CHKERRQ(ierr); 732 ierr = VecGetArrayRead(Xnet,&xnet);CHKERRQ(ierr); 733 734 /* Generator subsystem */ 735 for (i=0; i < ngen; i++) { 736 Eqp = xgen[idx]; 737 Edp = xgen[idx+1]; 738 delta = xgen[idx+2]; 739 Id = xgen[idx+4]; 740 Iq = xgen[idx+5]; 741 Efd = xgen[idx+6]; 742 743 /* fgen[idx] = (-Eqp - (Xd[i] - Xdp[i])*Id + Efd)/Td0p[i]; */ 744 row[0] = idx; 745 col[0] = idx; col[1] = idx+4; col[2] = idx+6; 746 val[0] = -1/ Td0p[i]; val[1] = -(Xd[i] - Xdp[i])/ Td0p[i]; val[2] = 1/Td0p[i]; 747 748 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 749 750 /* fgen[idx+1] = (-Edp + (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */ 751 row[0] = idx + 1; 752 col[0] = idx + 1; col[1] = idx+5; 753 val[0] = -1/Tq0p[i]; val[1] = (Xq[i] - Xqp[i])/Tq0p[i]; 754 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 755 756 /* fgen[idx+2] = w - w_s; */ 757 row[0] = idx + 2; 758 col[0] = idx + 2; col[1] = idx + 3; 759 val[0] = 0; val[1] = 1; 760 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 761 762 /* fgen[idx+3] = (TM[i] - Edp*Id - Eqp*Iq - (Xqp[i] - Xdp[i])*Id*Iq - D[i]*(w - w_s))/M[i]; */ 763 row[0] = idx + 3; 764 col[0] = idx; col[1] = idx + 1; col[2] = idx + 3; col[3] = idx + 4; col[4] = idx + 5; 765 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]; 766 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 767 768 Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */ 769 Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */ 770 ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr); 771 772 det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i]; 773 774 Zdq_inv[0] = Rs[i]/det; 775 Zdq_inv[1] = Xqp[i]/det; 776 Zdq_inv[2] = -Xdp[i]/det; 777 Zdq_inv[3] = Rs[i]/det; 778 779 dVd_dVr = PetscSinScalar(delta); dVd_dVi = -PetscCosScalar(delta); 780 dVq_dVr = PetscCosScalar(delta); dVq_dVi = PetscSinScalar(delta); 781 dVd_ddelta = Vr*PetscCosScalar(delta) + Vi*PetscSinScalar(delta); 782 dVq_ddelta = -Vr*PetscSinScalar(delta) + Vi*PetscCosScalar(delta); 783 784 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */ 785 row[0] = idx+4; 786 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 787 val[0] = -Zdq_inv[1]; val[1] = -Zdq_inv[0]; val[2] = Zdq_inv[0]*dVd_ddelta + Zdq_inv[1]*dVq_ddelta; 788 col[3] = idx + 4; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 789 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; 790 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 791 792 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */ 793 row[0] = idx+5; 794 col[0] = idx; col[1] = idx+1; col[2] = idx + 2; 795 val[0] = -Zdq_inv[3]; val[1] = -Zdq_inv[2]; val[2] = Zdq_inv[2]*dVd_ddelta + Zdq_inv[3]*dVq_ddelta; 796 col[3] = idx + 5; col[4] = net_start+2*gbus[i]; col[5] = net_start + 2*gbus[i]+1; 797 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; 798 ierr = MatSetValues(J,1,row,6,col,val,INSERT_VALUES);CHKERRQ(ierr); 799 800 dIGr_ddelta = Id*PetscCosScalar(delta) - Iq*PetscSinScalar(delta); 801 dIGi_ddelta = Id*PetscSinScalar(delta) + Iq*PetscCosScalar(delta); 802 dIGr_dId = PetscSinScalar(delta); dIGr_dIq = PetscCosScalar(delta); 803 dIGi_dId = -PetscCosScalar(delta); dIGi_dIq = PetscSinScalar(delta); 804 805 /* fnet[2*gbus[i]] -= IGi; */ 806 row[0] = net_start + 2*gbus[i]; 807 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 808 val[0] = -dIGi_ddelta; val[1] = -dIGi_dId; val[2] = -dIGi_dIq; 809 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 810 811 /* fnet[2*gbus[i]+1] -= IGr; */ 812 row[0] = net_start + 2*gbus[i]+1; 813 col[0] = idx+2; col[1] = idx + 4; col[2] = idx + 5; 814 val[0] = -dIGr_ddelta; val[1] = -dIGr_dId; val[2] = -dIGr_dIq; 815 ierr = MatSetValues(J,1,row,3,col,val,INSERT_VALUES);CHKERRQ(ierr); 816 817 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 818 819 /* fgen[idx+6] = (-KE[i]*Efd - SE + VR)/TE[i]; */ 820 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */ 821 822 dSE_dEfd = k1[i]*k2[i]*PetscExpScalar(k2[i]*Efd); 823 824 row[0] = idx + 6; 825 col[0] = idx + 6; col[1] = idx + 8; 826 val[0] = (-KE[i] - dSE_dEfd)/TE[i]; val[1] = 1/TE[i]; 827 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 828 829 /* Exciter differential equations */ 830 831 /* fgen[idx+7] = (-RF + KF[i]*Efd/TF[i])/TF[i]; */ 832 row[0] = idx + 7; 833 col[0] = idx + 6; col[1] = idx + 7; 834 val[0] = (KF[i]/TF[i])/TF[i]; val[1] = -1/TF[i]; 835 ierr = MatSetValues(J,1,row,2,col,val,INSERT_VALUES);CHKERRQ(ierr); 836 837 /* fgen[idx+8] = (-VR + KA[i]*RF - KA[i]*KF[i]*Efd/TF[i] + KA[i]*(Vref[i] - Vm))/TA[i]; */ 838 /* Vm = (Vd^2 + Vq^2)^0.5; */ 839 840 row[0] = idx + 8; 841 if (VRatmax[i]) { 842 col[0] = idx + 8; val[0] = 1.0; 843 ierr = MatSetValues(J,1,row,1,col,val,INSERT_VALUES);CHKERRQ(ierr); 844 } else if (VRatmin[i]) { 845 col[0] = idx + 8; val[0] = -1.0; 846 ierr = MatSetValues(J,1,row,1,col,val,INSERT_VALUES);CHKERRQ(ierr); 847 } else { 848 dVm_dVd = Vd/Vm; dVm_dVq = Vq/Vm; 849 dVm_dVr = dVm_dVd*dVd_dVr + dVm_dVq*dVq_dVr; 850 dVm_dVi = dVm_dVd*dVd_dVi + dVm_dVq*dVq_dVi; 851 row[0] = idx + 8; 852 col[0] = idx + 6; col[1] = idx + 7; col[2] = idx + 8; 853 val[0] = -(KA[i]*KF[i]/TF[i])/TA[i]; val[1] = KA[i]/TA[i]; val[2] = -1/TA[i]; 854 col[3] = net_start + 2*gbus[i]; col[4] = net_start + 2*gbus[i]+1; 855 val[3] = -KA[i]*dVm_dVr/TA[i]; val[4] = -KA[i]*dVm_dVi/TA[i]; 856 ierr = MatSetValues(J,1,row,5,col,val,INSERT_VALUES);CHKERRQ(ierr); 857 } 858 idx = idx + 9; 859 } 860 861 for (i=0; i<nbus; i++) { 862 ierr = MatGetRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 863 row[0] = net_start + 2*i; 864 for (k=0; k<ncols; k++) { 865 col[k] = net_start + cols[k]; 866 val[k] = yvals[k]; 867 } 868 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 869 ierr = MatRestoreRow(user->Ybus,2*i,&ncols,&cols,&yvals);CHKERRQ(ierr); 870 871 ierr = MatGetRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 872 row[0] = net_start + 2*i+1; 873 for (k=0; k<ncols; k++) { 874 col[k] = net_start + cols[k]; 875 val[k] = yvals[k]; 876 } 877 ierr = MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);CHKERRQ(ierr); 878 ierr = MatRestoreRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);CHKERRQ(ierr); 879 } 880 881 ierr = MatAssemblyBegin(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 882 ierr = MatAssemblyEnd(J,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 883 884 ierr = VecGetArrayRead(user->V0,&v0);CHKERRQ(ierr); 885 for (i=0; i < nload; i++) { 886 Vr = xnet[2*lbus[i]]; /* Real part of load bus voltage */ 887 Vi = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */ 888 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; Vm4 = Vm2*Vm2; 889 Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]); 890 PD = QD = 0.0; 891 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0; 892 for (k=0; k < ld_nsegsp[i]; k++) { 893 PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]); 894 dPD_dVr += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vr*PetscPowScalar(Vm,(ld_betap[k]-2)); 895 dPD_dVi += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vi*PetscPowScalar(Vm,(ld_betap[k]-2)); 896 } 897 for (k=0; k < ld_nsegsq[i]; k++) { 898 QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 899 dQD_dVr += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vr*PetscPowScalar(Vm,(ld_betaq[k]-2)); 900 dQD_dVi += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vi*PetscPowScalar(Vm,(ld_betaq[k]-2)); 901 } 902 903 /* IDr = (PD*Vr + QD*Vi)/Vm2; */ 904 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */ 905 906 dIDr_dVr = (dPD_dVr*Vr + dQD_dVr*Vi + PD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vr)/Vm4; 907 dIDr_dVi = (dPD_dVi*Vr + dQD_dVi*Vi + QD)/Vm2 - ((PD*Vr + QD*Vi)*2*Vi)/Vm4; 908 909 dIDi_dVr = (-dQD_dVr*Vr + dPD_dVr*Vi - QD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vr)/Vm4; 910 dIDi_dVi = (-dQD_dVi*Vr + dPD_dVi*Vi + PD)/Vm2 - ((-QD*Vr + PD*Vi)*2*Vi)/Vm4; 911 912 913 /* fnet[2*lbus[i]] += IDi; */ 914 row[0] = net_start + 2*lbus[i]; 915 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 916 val[0] = dIDi_dVr; val[1] = dIDi_dVi; 917 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 918 /* fnet[2*lbus[i]+1] += IDr; */ 919 row[0] = net_start + 2*lbus[i]+1; 920 col[0] = net_start + 2*lbus[i]; col[1] = net_start + 2*lbus[i]+1; 921 val[0] = dIDr_dVr; val[1] = dIDr_dVi; 922 ierr = MatSetValues(J,1,row,2,col,val,ADD_VALUES);CHKERRQ(ierr); 923 } 924 ierr = VecRestoreArrayRead(user->V0,&v0);CHKERRQ(ierr); 925 926 ierr = VecRestoreArrayRead(Xgen,&xgen);CHKERRQ(ierr); 927 ierr = VecRestoreArrayRead(Xnet,&xnet);CHKERRQ(ierr); 928 929 ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr); 930 931 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 932 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 933 PetscFunctionReturn(0); 934 } 935 936 /* 937 J = [I, 0 938 dg_dx, dg_dy] 939 */ 940 PetscErrorCode AlgJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx) 941 { 942 PetscErrorCode ierr; 943 Userctx *user=(Userctx*)ctx; 944 945 PetscFunctionBegin; 946 ierr = ResidualJacobian(X,A,B,ctx);CHKERRQ(ierr); 947 ierr = MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);CHKERRQ(ierr); 948 ierr = MatZeroRowsIS(A,user->is_diff,1.0,NULL,NULL);CHKERRQ(ierr); 949 PetscFunctionReturn(0); 950 } 951 952 /* 953 J = [-df_dx, -df_dy 954 dg_dx, dg_dy] 955 */ 956 957 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec X,Mat A,Mat B,void *ctx) 958 { 959 PetscErrorCode ierr; 960 Userctx *user=(Userctx*)ctx; 961 962 PetscFunctionBegin; 963 user->t = t; 964 965 ierr = ResidualJacobian(X,A,B,user);CHKERRQ(ierr); 966 967 PetscFunctionReturn(0); 968 } 969 970 /* 971 J = [df_dx-aI, df_dy 972 dg_dx, dg_dy] 973 */ 974 975 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx *user) 976 { 977 PetscErrorCode ierr; 978 PetscScalar atmp = (PetscScalar) a; 979 PetscInt i,row; 980 981 PetscFunctionBegin; 982 user->t = t; 983 984 ierr = RHSJacobian(ts,t,X,A,B,user);CHKERRQ(ierr); 985 ierr = MatScale(B,-1.0);CHKERRQ(ierr); 986 for (i=0;i < ngen;i++) { 987 row = 9*i; 988 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 989 row = 9*i+1; 990 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 991 row = 9*i+2; 992 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 993 row = 9*i+3; 994 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 995 row = 9*i+6; 996 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 997 row = 9*i+7; 998 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 999 row = 9*i+8; 1000 ierr = MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);CHKERRQ(ierr); 1001 } 1002 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1003 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1004 PetscFunctionReturn(0); 1005 } 1006 1007 int main(int argc,char **argv) 1008 { 1009 TS ts; 1010 SNES snes_alg; 1011 PetscErrorCode ierr; 1012 PetscMPIInt size; 1013 Userctx user; 1014 PetscViewer Xview,Ybusview,viewer; 1015 Vec X,F_alg; 1016 Mat J,A; 1017 PetscInt i,idx,*idx2; 1018 Vec Xdot; 1019 PetscScalar *x,*mat,*amat; 1020 const PetscScalar *rmat; 1021 Vec vatol; 1022 PetscInt *direction; 1023 PetscBool *terminate; 1024 const PetscInt *idx3; 1025 PetscScalar *vatoli; 1026 PetscInt k; 1027 1028 1029 ierr = PetscInitialize(&argc,&argv,"petscoptions",help);if (ierr) return ierr; 1030 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 1031 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 1032 1033 user.neqs_gen = 9*ngen; /* # eqs. for generator subsystem */ 1034 user.neqs_net = 2*nbus; /* # eqs. for network subsystem */ 1035 user.neqs_pgrid = user.neqs_gen + user.neqs_net; 1036 1037 /* Create indices for differential and algebraic equations */ 1038 1039 ierr = PetscMalloc1(7*ngen,&idx2);CHKERRQ(ierr); 1040 for (i=0; i<ngen; i++) { 1041 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; 1042 idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8; 1043 } 1044 ierr = ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff);CHKERRQ(ierr); 1045 ierr = ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg);CHKERRQ(ierr); 1046 ierr = PetscFree(idx2);CHKERRQ(ierr); 1047 1048 /* Read initial voltage vector and Ybus */ 1049 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"X.bin",FILE_MODE_READ,&Xview);CHKERRQ(ierr); 1050 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ybus.bin",FILE_MODE_READ,&Ybusview);CHKERRQ(ierr); 1051 1052 ierr = VecCreate(PETSC_COMM_WORLD,&user.V0);CHKERRQ(ierr); 1053 ierr = VecSetSizes(user.V0,PETSC_DECIDE,user.neqs_net);CHKERRQ(ierr); 1054 ierr = VecLoad(user.V0,Xview);CHKERRQ(ierr); 1055 1056 ierr = MatCreate(PETSC_COMM_WORLD,&user.Ybus);CHKERRQ(ierr); 1057 ierr = MatSetSizes(user.Ybus,PETSC_DECIDE,PETSC_DECIDE,user.neqs_net,user.neqs_net);CHKERRQ(ierr); 1058 ierr = MatSetType(user.Ybus,MATBAIJ);CHKERRQ(ierr); 1059 /* ierr = MatSetBlockSize(user.Ybus,2);CHKERRQ(ierr); */ 1060 ierr = MatLoad(user.Ybus,Ybusview);CHKERRQ(ierr); 1061 1062 /* Set run time options */ 1063 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");CHKERRQ(ierr); 1064 { 1065 user.tfaulton = 1.0; 1066 user.tfaultoff = 1.2; 1067 user.Rfault = 0.0001; 1068 user.setisdiff = PETSC_FALSE; 1069 user.semiexplicit = PETSC_FALSE; 1070 user.faultbus = 8; 1071 ierr = PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL);CHKERRQ(ierr); 1072 ierr = PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL);CHKERRQ(ierr); 1073 ierr = PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL);CHKERRQ(ierr); 1074 user.t0 = 0.0; 1075 user.tmax = 5.0; 1076 ierr = PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL);CHKERRQ(ierr); 1077 ierr = PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL);CHKERRQ(ierr); 1078 ierr = PetscOptionsBool("-setisdiff","","",user.setisdiff,&user.setisdiff,NULL);CHKERRQ(ierr); 1079 ierr = PetscOptionsBool("-dae_semiexplicit","","",user.semiexplicit,&user.semiexplicit,NULL);CHKERRQ(ierr); 1080 } 1081 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1082 1083 ierr = PetscViewerDestroy(&Xview);CHKERRQ(ierr); 1084 ierr = PetscViewerDestroy(&Ybusview);CHKERRQ(ierr); 1085 1086 /* Create DMs for generator and network subsystems */ 1087 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); 1088 ierr = DMSetOptionsPrefix(user.dmgen,"dmgen_");CHKERRQ(ierr); 1089 ierr = DMSetFromOptions(user.dmgen);CHKERRQ(ierr); 1090 ierr = DMSetUp(user.dmgen);CHKERRQ(ierr); 1091 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); 1092 ierr = DMSetOptionsPrefix(user.dmnet,"dmnet_");CHKERRQ(ierr); 1093 ierr = DMSetFromOptions(user.dmnet);CHKERRQ(ierr); 1094 ierr = DMSetUp(user.dmnet);CHKERRQ(ierr); 1095 /* Create a composite DM packer and add the two DMs */ 1096 ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid);CHKERRQ(ierr); 1097 ierr = DMSetOptionsPrefix(user.dmpgrid,"pgrid_");CHKERRQ(ierr); 1098 ierr = DMCompositeAddDM(user.dmpgrid,user.dmgen);CHKERRQ(ierr); 1099 ierr = DMCompositeAddDM(user.dmpgrid,user.dmnet);CHKERRQ(ierr); 1100 1101 ierr = DMCreateGlobalVector(user.dmpgrid,&X);CHKERRQ(ierr); 1102 1103 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 1104 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,user.neqs_pgrid,user.neqs_pgrid);CHKERRQ(ierr); 1105 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 1106 ierr = PreallocateJacobian(J,&user);CHKERRQ(ierr); 1107 1108 /* Create matrix to save solutions at each time step */ 1109 user.stepnum = 0; 1110 1111 ierr = MatCreateSeqDense(PETSC_COMM_SELF,user.neqs_pgrid+1,1002,NULL,&user.Sol);CHKERRQ(ierr); 1112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1113 Create timestepping solver context 1114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1115 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 1116 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 1117 if (user.semiexplicit) { 1118 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr); 1119 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); 1120 ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,&user);CHKERRQ(ierr); 1121 } else { 1122 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 1123 ierr = TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);CHKERRQ(ierr); 1124 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&user);CHKERRQ(ierr); 1125 ierr = TSSetIJacobian(ts,J,J,(TSIJacobian)IJacobian,&user);CHKERRQ(ierr); 1126 } 1127 ierr = TSSetApplicationContext(ts,&user);CHKERRQ(ierr); 1128 1129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1130 Set initial conditions 1131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1132 ierr = SetInitialGuess(X,&user);CHKERRQ(ierr); 1133 /* Just to set up the Jacobian structure */ 1134 1135 ierr = VecDuplicate(X,&Xdot);CHKERRQ(ierr); 1136 ierr = IJacobian(ts,0.0,X,Xdot,0.0,J,J,&user);CHKERRQ(ierr); 1137 ierr = VecDestroy(&Xdot);CHKERRQ(ierr); 1138 1139 /* Save initial solution */ 1140 1141 idx=user.stepnum*(user.neqs_pgrid+1); 1142 ierr = MatDenseGetArray(user.Sol,&mat);CHKERRQ(ierr); 1143 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 1144 1145 mat[idx] = 0.0; 1146 1147 ierr = PetscArraycpy(mat+idx+1,x,user.neqs_pgrid);CHKERRQ(ierr); 1148 ierr = MatDenseRestoreArray(user.Sol,&mat);CHKERRQ(ierr); 1149 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 1150 user.stepnum++; 1151 1152 ierr = TSSetMaxTime(ts,user.tmax);CHKERRQ(ierr); 1153 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 1154 ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr); 1155 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 1156 ierr = TSSetPostStep(ts,SaveSolution);CHKERRQ(ierr); 1157 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 1158 1159 ierr = PetscMalloc1((2*ngen+2),&direction);CHKERRQ(ierr); 1160 ierr = PetscMalloc1((2*ngen+2),&terminate);CHKERRQ(ierr); 1161 direction[0] = direction[1] = 1; 1162 terminate[0] = terminate[1] = PETSC_FALSE; 1163 for (i=0; i < ngen;i++) { 1164 direction[2+2*i] = -1; direction[2+2*i+1] = 1; 1165 terminate[2+2*i] = terminate[2+2*i+1] = PETSC_FALSE; 1166 } 1167 1168 ierr = TSSetEventHandler(ts,2*ngen+2,direction,terminate,EventFunction,PostEventFunction,(void*)&user);CHKERRQ(ierr); 1169 1170 if (user.semiexplicit) { 1171 /* Use a semi-explicit approach with the time-stepping done by an explicit method and the 1172 algrebraic part solved via PostStage and PostEvaluate callbacks 1173 */ 1174 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr); 1175 ierr = TSSetPostStage(ts,PostStage);CHKERRQ(ierr); 1176 ierr = TSSetPostEvaluate(ts,PostEvaluate);CHKERRQ(ierr); 1177 } 1178 1179 1180 if (user.setisdiff) { 1181 /* Create vector of absolute tolerances and set the algebraic part to infinity */ 1182 ierr = VecDuplicate(X,&vatol);CHKERRQ(ierr); 1183 ierr = VecSet(vatol,100000.0);CHKERRQ(ierr); 1184 ierr = VecGetArray(vatol,&vatoli);CHKERRQ(ierr); 1185 ierr = ISGetIndices(user.is_diff,&idx3);CHKERRQ(ierr); 1186 for (k=0; k < 7*ngen; k++) vatoli[idx3[k]] = 1e-2; 1187 ierr = VecRestoreArray(vatol,&vatoli);CHKERRQ(ierr); 1188 } 1189 1190 /* Create the nonlinear solver for solving the algebraic system */ 1191 /* Note that although the algebraic system needs to be solved only for 1192 Idq and V, we reuse the entire system including xgen. The xgen 1193 variables are held constant by setting their residuals to 0 and 1194 putting a 1 on the Jacobian diagonal for xgen rows 1195 */ 1196 1197 ierr = VecDuplicate(X,&F_alg);CHKERRQ(ierr); 1198 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_alg);CHKERRQ(ierr); 1199 ierr = SNESSetFunction(snes_alg,F_alg,AlgFunction,&user);CHKERRQ(ierr); 1200 ierr = SNESSetJacobian(snes_alg,J,J,AlgJacobian,&user);CHKERRQ(ierr); 1201 ierr = SNESSetFromOptions(snes_alg);CHKERRQ(ierr); 1202 1203 user.snes_alg=snes_alg; 1204 1205 /* Solve */ 1206 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1207 1208 ierr = MatAssemblyBegin(user.Sol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1209 ierr = MatAssemblyEnd(user.Sol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1210 1211 ierr = MatCreateSeqDense(PETSC_COMM_SELF,user.neqs_pgrid+1,user.stepnum,NULL,&A);CHKERRQ(ierr); 1212 ierr = MatDenseGetArrayRead(user.Sol,&rmat);CHKERRQ(ierr); 1213 ierr = MatDenseGetArray(A,&amat);CHKERRQ(ierr); 1214 ierr = PetscArraycpy(amat,rmat,user.stepnum*(user.neqs_pgrid+1));CHKERRQ(ierr); 1215 ierr = MatDenseRestoreArray(A,&amat);CHKERRQ(ierr); 1216 ierr = MatDenseRestoreArrayRead(user.Sol,&rmat);CHKERRQ(ierr); 1217 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"out.bin",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 1218 ierr = MatView(A,viewer);CHKERRQ(ierr); 1219 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1220 ierr = MatDestroy(&A);CHKERRQ(ierr); 1221 1222 ierr = PetscFree(direction);CHKERRQ(ierr); 1223 ierr = PetscFree(terminate);CHKERRQ(ierr); 1224 ierr = SNESDestroy(&snes_alg);CHKERRQ(ierr); 1225 ierr = VecDestroy(&F_alg);CHKERRQ(ierr); 1226 ierr = MatDestroy(&J);CHKERRQ(ierr); 1227 ierr = MatDestroy(&user.Ybus);CHKERRQ(ierr); 1228 ierr = MatDestroy(&user.Sol);CHKERRQ(ierr); 1229 ierr = VecDestroy(&X);CHKERRQ(ierr); 1230 ierr = VecDestroy(&user.V0);CHKERRQ(ierr); 1231 ierr = DMDestroy(&user.dmgen);CHKERRQ(ierr); 1232 ierr = DMDestroy(&user.dmnet);CHKERRQ(ierr); 1233 ierr = DMDestroy(&user.dmpgrid);CHKERRQ(ierr); 1234 ierr = ISDestroy(&user.is_diff);CHKERRQ(ierr); 1235 ierr = ISDestroy(&user.is_alg);CHKERRQ(ierr); 1236 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1237 if (user.setisdiff) { 1238 ierr = VecDestroy(&vatol);CHKERRQ(ierr); 1239 } 1240 ierr = PetscFinalize(); 1241 return ierr; 1242 } 1243 1244 /*TEST 1245 1246 build: 1247 requires: double !complex !define(PETSC_USE_64BIT_INDICES) 1248 1249 test: 1250 suffix: implicit 1251 args: -ts_monitor -snes_monitor_short 1252 localrunfiles: petscoptions X.bin Ybus.bin 1253 1254 test: 1255 suffix: semiexplicit 1256 args: -ts_monitor -snes_monitor_short -dae_semiexplicit -ts_rk_type 2a 1257 localrunfiles: petscoptions X.bin Ybus.bin 1258 1259 test: 1260 suffix: steprestart 1261 # needs ARKIMEX methods with all implicit stages since the mass matrix is not the identity 1262 args: -ts_monitor -snes_monitor_short -ts_type arkimex -ts_arkimex_type prssp2 1263 localrunfiles: petscoptions X.bin Ybus.bin 1264 1265 TEST*/ 1266