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