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