1 2 static char help[] = "Adjoint and tangent linear sensitivity analysis of the basic equation for generator stability analysis.\n"; 3 4 /*F 5 6 \begin{eqnarray} 7 \frac{d \theta}{dt} = \omega_b (\omega - \omega_s) 8 \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\ 9 \end{eqnarray} 10 11 F*/ 12 13 /* 14 This code demonstrate the sensitivity analysis interface to a system of ordinary differential equations with discontinuities. 15 It computes the sensitivities of an integral cost function 16 \int c*max(0,\theta(t)-u_s)^beta dt 17 w.r.t. initial conditions and the parameter P_m. 18 Backward Euler method is used for time integration. 19 The discontinuities are detected with TSEvent. 20 */ 21 22 #include <petscts.h> 23 #include "ex3.h" 24 25 int main(int argc,char **argv) 26 { 27 TS ts,quadts; /* ODE integrator */ 28 Vec U; /* solution will be stored here */ 29 PetscMPIInt size; 30 PetscInt n = 2; 31 AppCtx ctx; 32 PetscScalar *u; 33 PetscReal du[2] = {0.0,0.0}; 34 PetscBool ensemble = PETSC_FALSE,flg1,flg2; 35 PetscReal ftime; 36 PetscInt steps; 37 PetscScalar *x_ptr,*y_ptr,*s_ptr; 38 Vec lambda[1],q,mu[1]; 39 PetscInt direction[2]; 40 PetscBool terminate[2]; 41 Mat qgrad; 42 Mat sp; /* Forward sensitivity matrix */ 43 SAMethod sa; 44 45 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 46 Initialize program 47 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 48 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 49 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 50 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create necessary matrix and vectors 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 PetscCall(MatCreate(PETSC_COMM_WORLD,&ctx.Jac)); 56 PetscCall(MatSetSizes(ctx.Jac,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 57 PetscCall(MatSetType(ctx.Jac,MATDENSE)); 58 PetscCall(MatSetFromOptions(ctx.Jac)); 59 PetscCall(MatSetUp(ctx.Jac)); 60 PetscCall(MatCreateVecs(ctx.Jac,&U,NULL)); 61 PetscCall(MatCreate(PETSC_COMM_WORLD,&ctx.Jacp)); 62 PetscCall(MatSetSizes(ctx.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1)); 63 PetscCall(MatSetFromOptions(ctx.Jacp)); 64 PetscCall(MatSetUp(ctx.Jacp)); 65 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&ctx.DRDP)); 66 PetscCall(MatSetUp(ctx.DRDP)); 67 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&ctx.DRDU)); 68 PetscCall(MatSetUp(ctx.DRDU)); 69 70 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 71 Set runtime options 72 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 73 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options",""); 74 { 75 ctx.beta = 2; 76 ctx.c = 10000.0; 77 ctx.u_s = 1.0; 78 ctx.omega_s = 1.0; 79 ctx.omega_b = 120.0*PETSC_PI; 80 ctx.H = 5.0; 81 PetscCall(PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL)); 82 ctx.D = 5.0; 83 PetscCall(PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL)); 84 ctx.E = 1.1378; 85 ctx.V = 1.0; 86 ctx.X = 0.545; 87 ctx.Pmax = ctx.E*ctx.V/ctx.X; 88 ctx.Pmax_ini = ctx.Pmax; 89 PetscCall(PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL)); 90 ctx.Pm = 1.1; 91 PetscCall(PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL)); 92 ctx.tf = 0.1; 93 ctx.tcl = 0.2; 94 PetscCall(PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL)); 95 PetscCall(PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL)); 96 PetscCall(PetscOptionsBool("-ensemble","Run ensemble of different initial conditions","",ensemble,&ensemble,NULL)); 97 if (ensemble) { 98 ctx.tf = -1; 99 ctx.tcl = -1; 100 } 101 102 PetscCall(VecGetArray(U,&u)); 103 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax); 104 u[1] = 1.0; 105 PetscCall(PetscOptionsRealArray("-u","Initial solution","",u,&n,&flg1)); 106 n = 2; 107 PetscCall(PetscOptionsRealArray("-du","Perturbation in initial solution","",du,&n,&flg2)); 108 u[0] += du[0]; 109 u[1] += du[1]; 110 PetscCall(VecRestoreArray(U,&u)); 111 if (flg1 || flg2) { 112 ctx.tf = -1; 113 ctx.tcl = -1; 114 } 115 sa = SA_ADJ; 116 PetscCall(PetscOptionsEnum("-sa_method","Sensitivity analysis method (adj or tlm)","",SAMethods,(PetscEnum)sa,(PetscEnum*)&sa,NULL)); 117 } 118 PetscOptionsEnd(); 119 120 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 121 Create timestepping solver context 122 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 123 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 124 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 125 PetscCall(TSSetType(ts,TSBEULER)); 126 PetscCall(TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx)); 127 PetscCall(TSSetRHSJacobian(ts,ctx.Jac,ctx.Jac,(TSRHSJacobian)RHSJacobian,&ctx)); 128 129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 Set initial conditions 131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 132 PetscCall(TSSetSolution(ts,U)); 133 134 /* Set RHS JacobianP */ 135 PetscCall(TSSetRHSJacobianP(ts,ctx.Jacp,RHSJacobianP,&ctx)); 136 137 PetscCall(TSCreateQuadratureTS(ts,PETSC_FALSE,&quadts)); 138 PetscCall(TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx)); 139 PetscCall(TSSetRHSJacobian(quadts,ctx.DRDU,ctx.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx)); 140 PetscCall(TSSetRHSJacobianP(quadts,ctx.DRDP,DRDPJacobianTranspose,&ctx)); 141 if (sa == SA_ADJ) { 142 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 Save trajectory of solution so that TSAdjointSolve() may be used 144 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 145 PetscCall(TSSetSaveTrajectory(ts)); 146 PetscCall(MatCreateVecs(ctx.Jac,&lambda[0],NULL)); 147 PetscCall(MatCreateVecs(ctx.Jacp,&mu[0],NULL)); 148 PetscCall(TSSetCostGradients(ts,1,lambda,mu)); 149 } 150 151 if (sa == SA_TLM) { 152 PetscScalar val[2]; 153 PetscInt row[]={0,1},col[]={0}; 154 155 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&qgrad)); 156 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&sp)); 157 PetscCall(TSForwardSetSensitivities(ts,1,sp)); 158 PetscCall(TSForwardSetSensitivities(quadts,1,qgrad)); 159 val[0] = 1./PetscSqrtScalar(1.-(ctx.Pm/ctx.Pmax)*(ctx.Pm/ctx.Pmax))/ctx.Pmax; 160 val[1] = 0.0; 161 PetscCall(MatSetValues(sp,2,row,1,col,val,INSERT_VALUES)); 162 PetscCall(MatAssemblyBegin(sp,MAT_FINAL_ASSEMBLY)); 163 PetscCall(MatAssemblyEnd(sp,MAT_FINAL_ASSEMBLY)); 164 } 165 166 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 167 Set solver options 168 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 169 PetscCall(TSSetMaxTime(ts,1.0)); 170 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 171 PetscCall(TSSetTimeStep(ts,0.03125)); 172 PetscCall(TSSetFromOptions(ts)); 173 174 direction[0] = direction[1] = 1; 175 terminate[0] = terminate[1] = PETSC_FALSE; 176 177 PetscCall(TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&ctx)); 178 179 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 180 Solve nonlinear system 181 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 182 if (ensemble) { 183 for (du[1] = -2.5; du[1] <= .01; du[1] += .1) { 184 PetscCall(VecGetArray(U,&u)); 185 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax); 186 u[1] = ctx.omega_s; 187 u[0] += du[0]; 188 u[1] += du[1]; 189 PetscCall(VecRestoreArray(U,&u)); 190 PetscCall(TSSetTimeStep(ts,0.03125)); 191 PetscCall(TSSolve(ts,U)); 192 } 193 } else { 194 PetscCall(TSSolve(ts,U)); 195 } 196 PetscCall(TSGetSolveTime(ts,&ftime)); 197 PetscCall(TSGetStepNumber(ts,&steps)); 198 199 if (sa == SA_ADJ) { 200 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 201 Adjoint model starts here 202 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 203 /* Set initial conditions for the adjoint integration */ 204 PetscCall(VecGetArray(lambda[0],&y_ptr)); 205 y_ptr[0] = 0.0; y_ptr[1] = 0.0; 206 PetscCall(VecRestoreArray(lambda[0],&y_ptr)); 207 208 PetscCall(VecGetArray(mu[0],&x_ptr)); 209 x_ptr[0] = 0.0; 210 PetscCall(VecRestoreArray(mu[0],&x_ptr)); 211 212 PetscCall(TSAdjointSolve(ts)); 213 214 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n lambda: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n")); 215 PetscCall(VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD)); 216 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n mu: d[Psi(tf)]/d[pm]\n")); 217 PetscCall(VecView(mu[0],PETSC_VIEWER_STDOUT_WORLD)); 218 PetscCall(TSGetCostIntegral(ts,&q)); 219 PetscCall(VecGetArray(q,&x_ptr)); 220 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(x_ptr[0]-ctx.Pm))); 221 PetscCall(VecRestoreArray(q,&x_ptr)); 222 PetscCall(ComputeSensiP(lambda[0],mu[0],&ctx)); 223 PetscCall(VecGetArray(mu[0],&x_ptr)); 224 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)x_ptr[0])); 225 PetscCall(VecRestoreArray(mu[0],&x_ptr)); 226 PetscCall(VecDestroy(&lambda[0])); 227 PetscCall(VecDestroy(&mu[0])); 228 } 229 if (sa == SA_TLM) { 230 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n trajectory sensitivity: d[phi(tf)]/d[pm] d[omega(tf)]/d[pm]\n")); 231 PetscCall(MatView(sp,PETSC_VIEWER_STDOUT_WORLD)); 232 PetscCall(TSGetCostIntegral(ts,&q)); 233 PetscCall(VecGetArray(q,&s_ptr)); 234 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(s_ptr[0]-ctx.Pm))); 235 PetscCall(VecRestoreArray(q,&s_ptr)); 236 PetscCall(MatDenseGetArray(qgrad,&s_ptr)); 237 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)s_ptr[0])); 238 PetscCall(MatDenseRestoreArray(qgrad,&s_ptr)); 239 PetscCall(MatDestroy(&qgrad)); 240 PetscCall(MatDestroy(&sp)); 241 } 242 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243 Free work space. All PETSc objects should be destroyed when they are no longer needed. 244 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 245 PetscCall(MatDestroy(&ctx.Jac)); 246 PetscCall(MatDestroy(&ctx.Jacp)); 247 PetscCall(MatDestroy(&ctx.DRDU)); 248 PetscCall(MatDestroy(&ctx.DRDP)); 249 PetscCall(VecDestroy(&U)); 250 PetscCall(TSDestroy(&ts)); 251 PetscCall(PetscFinalize()); 252 return 0; 253 } 254 255 /*TEST 256 257 build: 258 requires: !complex !single 259 260 test: 261 args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu 262 263 test: 264 suffix: 2 265 args: -sa_method tlm -ts_type cn -pc_type lu 266 267 test: 268 suffix: 3 269 args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp 270 271 TEST*/ 272